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Lift/Drag Ratios of Optimised Slewed Elliptic 
Wings at Supersonic Speeds 


J. H. B. SMITH, M.Sc. 


(Royal Aircraft Establishment, Farnborough) 


SUMMARY: Previous work by R. T. Jones on the drag minimisation of elliptic 

wings is extended to the case of the slewed wing with thickness. These results 

are used to calculate lift/drag ratios of idealised configurations related to a 

supersonic transport aircraft. The values of lift/drag ratio and optimum 

slenderness ratio found are comparable with those calculated earlier in studies 
of delta-like plan forms. 


1. Introduction 


This paper gives the results of calculations of the lift/drag ratio for wings of 
elliptic plan form flying at supersonic speeds with their axes of symmetry slewed* 
with respect to the flight direction. It is assumed that the disturbances produced 
by the wing are small, that the Mach number is neither very large nor very near unity 
and that separation takes place from the trailing edge (A BC in Fig. 1) of the wing 
only, the wake remaining sensibly flat so long as it affects the wing. Under these 
assumptions linearised supersonic thin-wing theory can be used and volume and 
lifting effects are then separable. The volume of the wing is distributed over the 
plan form in such a way that the drag due to volume is minimised according to 
thin-wing theory. (The thickness is assumed to vanish at the edges.) The wing is 
warped so that at a design lift coefficient the lift-dependent drag is a minimum also. 


We have not considered whether the postulated type of flow is likely to be 
achieved in a real fluid; nor have we attempted to show that the configurations are 
potentially usable for aircraft. The calculations are merely intended to bring out 
the effects of certain geometrical parameters which are aerodynamically significant, 
and to be generally comparable with those reported by Kichemann*’ for other 
plan forms. 


*The term “slew” is used to denote a rotation about the yaw axis which transforms a basic 
shape into the plan form under consideration. The basic shape is here an ellipse with its 
major axis in the flight direction, so that the aspect ratio is a minimum at zero slew. This 
seems reasonable in the present discussion of supersonic performance; but it can be main- 
tained that a configuration with zero slew should be of high aspect ratio, like a configuration 
with zero sweepback. The term “ yaw” is reserved for rotations of the resulting wing away 
from it direction of steady flight 
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The work was done for two reasons. Various systems of variable geometry 
have been proposed to avoid the compromise in the design of a supersonic aircraft 
between the span-to-length ratios required for efficient supersonic flight and those 
for conventional take-off and landing. An all-wing aircraft which could fly with its 
axis at different angles of slew to the flight direction would embody a system of 
variable geometry of great mechanical simplicity. A plan form with a continuously 
curved edge is an obvious advantage for such an aircraft and an ellipse offers 
sufficient degrees of freedom for a preliminary investigation. One purpose of the 
present work was to discover what angles of slew would be desirable for the cruising 
flight of a supersonic transport aircraft of elliptic plan form and how the level of 
lift/drag ratio which might be attained would compare with that of a laterally 
symmetric aircraft. There has also been considerable interest in recent years in the 
span-to-length ratios suitable for delta-like aircraft designed for efficient supersonic 
flight. A complete investigation of the aerodynamics of such aircraft is only feasible 
if slender thin-wing theory is used (i.e. the limiting form of the theory of this 
paper as A (M*-1)' * tends to zero). The present calculations for slewed ellipses 
provide evidence which is relevant in a general sense to the choice of span-to-length 
ratio for a delta-like plan form. 


The choice of the ellipse as the plan form for these calculations is dictated by 
its mathematical simplicity. This makes it possible to express the minimum drag 
of the wing with given lift or volume in a simple closed form. For the slewed wing 
of given lift and zero volume and for the laterally symmetric uncambered wing of 
given volume, the expressions have been given by Jones”. The extension to the 
slewed wing of given volume is given in the Appendix to this paper. 


Lift/drag ratios for slewed ellipses have been quoted in a paper? by Jones. 
It is not clear from this paper whether the drag due to volume has been minimised 
and a fuller investigation was felt to be necessary. 


NOTATION 
A aspect ratio 
a,b  semi-minor axis and semi-major axis. respectively 
a see Fig. | and below equation (3) 
b’ the semi-span of the yawed ellipse (see Fig. | and below equation (3)) 
C,, drag coefficient 
friction drag coefficient 
C,, drag coefficient of uncambered wing at zero lift 
Cy,« Wave drag coefficient of uncambered wing at zero lift 
C, lift coefficient 
D drag 
D(#) drag of oblique area distribution (see Appendix) 
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d_ half length of oblique area distribution (see Appendix) 
K a lift-dependent drag factor 
K, vortex drag factor (equation (11)) 
lift-dependent wave-drag factor (equation (12)) 
K, zero-lift wave-drag factor (equation (10)) 
Lsiift 
/ overall length of wing 
M Mach number 
m= cos 4 (see Appendix) 
m’ see Fig. 1 and below equation (3) 
p plan form parameter (wing area/(span x length)) 
q_ kinetic pressure 
R_ Reynolds number 
A‘Z} real part of Z 
S wing area 
s half the overall span 
¢ maximum wing thickness 
U_ free-stream velocity 
VY volume 


v.¥,zZ right-handed rectangular axes, origin at centre of plan form, Ox 
streamwise, Oy to starboard 


B= /(M?*-1) 

see Appendix 

u Mach angle, =cosec~* M 

volume parameter, 


v~ angle of slew of major axis to flight direction 


Suffixes 


cr cruising conditions 


opt at C, for maximu.n L/D 


2. Assumptions and Formulae 


According to thin-wing theory the pressure drag of a wing of finite volume 
carrying a certain lift is the sum of the drag at zero lift of a wing having top and 
bottom symmetry with the same thickness distribution and of the drag of the 
camber surface of the wing at the same lift. It is approximately true that the 
skin-friction drag of a wing is independent of a small amount of warp applied to it 
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provided there is no separation ahead of the trailing edge. Hence, if we consider a 
wing warped so as to produce the least lift-dependent drag consistent with a given 
lift coefficient C,, its drag coefficient is 


= 


Here Cy, is the drag coefficient of an uncambered wing of the same volume 
distribution, 


where Cp,» is the zero-lift wave-drag coefficient and C,; is the coefficient of skin- 
friction drag. In equation (1), A is the aspect ratio and K is a lift-dependent drag 
factor for the camber surface alone. From Ref. 2, 


K a FW E |}. (3) 


for a slewed ellipse of minimum lift-dependent drag. 


Here £?=M?-1, b’ =(a? cos? 1+ b? sin? y)'/?, 


‘’=ab/b’, nv =(b? —a’) cos 


where 2a is the length of the minor axis of the ellipse, whose major axis of length 
2b is slewed at an angle ¥ to the stream of Mach number M. The geometry is 
shown in Fig. 1. The overall length of the plan form is 


(a? sin? b? cos? =2 + m’b’?)!/2, (4) 


From equation (23) of the Appendix, we find 


a a 
( m’ +i ( mn’ + 2i 
| 
| |e ( m | 


for the minimum wave-drag coefficient of a symmetric elliptic wing of given 
volume. The maximum thickness is ¢ and the volume is given by equation (20) as 


The volume parameter 7 is therefore given by 


where S is the plan form area. For the friction drag, flat-plate skin friction values 
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FiGuRE 1. Notation. 


for a turbulent boundary layer are taken and it is assumed that the friction acts on 
twice the plan area. The Reynolds number is based on the geometric mean chord 
and an altitude variation with Mach number is assumed such that the Reynolds 
number per foot is 3x 10°/M. If the overall length is 7 ft. and the plan form 
parameter (wing area)/(span x length) is p, the geometric mean chord is pl ft. and 
the Reynolds number is 3p/ = 10°/M. We use charts of skin-friction drag calcula- 
ted by the intermediate enthalpy method, arbitrarily choosing to take the wall 
enthalpy equal to the mean of the stream and recovery enthalpies. Simple curves 
are then fitted to the charts in the relevant range (1 = M=<5, 10° =R=<10"), 
this giving 


Cyt =(2:71 —0°333M) | 3-6 —0°597 log +0-0424 log 10-*, (8) 
i 


where the logarithms are natural and the mean skin friction coefficient has been 
assumed to be 1:25 times the local value. To fix the scale for the skin friction, a 
plan form area of 6,000 square feet is considered. 
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Non-aerodynamic considerations usually make it desirable to cruise at a lift 
coefficient below the aerodynamic optimum. A convenient value of cruising lift 
coefficient near those commonly found is 


=Crop:/ V2. 


The optimum lift coefficient is that for maximum lift/drag ratio and therefore the 


cruising lift/drag ratio is 
(L/D)x= *, 


Equations (2) to (9) determine the lift/drag ratio at cruise in terms of the 
volume and the plan form geometry. 


Of more general significance than these values of lift/drag ratio are the values 
of the drag factors K,, K, and K,, introduced in Ref. 1. These refer the volume 
wave-drag, the vortex drag and the lift-dependent wave-drag to certain standard 
values. K, is the ratio of the wave drag due to volume to that of a Sears-Haack 
body of the same length and volume. [The Sears-Haack body is an optimum 
according to slender-body theory in the sense that its drag is less than that of any 
other body of the same length / and volume V whose cross-sectional area S (x) 
satisfies § (0)=0=S ()=S’ (0)=S’ (I). Its drag, D, is given by 


128V? 
q 


K, is the ratio of the vortex drag to that of an elliptically loaded wing of the 
same span and lift. K,, is the ratio of the lift-dependent wave-drag to that found 
by not-so-slender wing theory“ for a wing with the same length and an elliptic 
distribution of cross-load. The lift-dependent drag factor, K, of equations (1) and 
(3) is given by 

K=K,+2(8s/l)’Ky. 


For the slewed elliptic wing considered here, we find 


| B?- +i ( 2i | 
(m +i?) | 
K,=1 (11) 
Ky= + a) Ea (mv +i |} -1| ‘ (12) 


Numerical calculations of lift/drag ratio at cruise and of K, and Ky have been 
made; the variation of these quantities is discussed in the following Section. 
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3. Discussion 


3.1. DraG Factors 


It is known from Ref. 2 that the lift-dependent drag of the slewed ellipse is a 
minimum when the lift is uniformly distributed over the plan form. The ellipse 
has the property that the lengths of any family of parallel chords are distributed 
elliptically. Hence the spanwise loading is elliptic and the vortex drag factor, Ky, 
is equal to one. 


The variation of the volume wave-drag factor K, and the lift-dependent wave- 
drag factor Ky is most easily explained in terms of the supersonic area rule (see 
Ref. 5, for example). This expresses the drag of a distribution of singularities, 
planar in the present case, as an average of the drags of a set of lineal singularity 
distributions. Each lineal distribution is obtained by considering a family of 
parallel Mach planes (i.e. planes inclined to the stream at the Mach angle, 
u=cosec~' M) and by transferring the singularities which lie in each plane to the 
point where it cuts a streamwise axis. Each family of Mach planes cuts the wing 
plane in a family of parallel lines, which make an angle with the stream direction 
lying between the Mach angle and a right angle. Thus the length of one of these 
lineal distributions is determined by the streamwise distance between two parallel 
tangents to the plan form. In Fig. 2(a) are shown the intersections of some members 
of one family of parallel planes with the wing plane. The singularities representing 
the volume and lift of the wing which are located along the chords P,Q, and P,Q, 
are transferred to R, and R, to form a lineal distribution between A and B. In 
Fig. 2(b) are shown the end-points A,, B, of some of the lineal distributions 
associated with a Mach number M=cosec ». A, B, is the longest and corresponds 
to the family of planes for which the lines make the Mach angle on one side of the 
stream direction. A, B, corresponds to the family of planes for which the lines are 
normal to the stream. A, B, corresponds to the planes for which the lines make 
the Mach angle on the other side of the stream direction; note that it is not 
necessarily the shortest. 


From Ref. 2, the lift distribution is uniform over the plan form when the 
lift-dependent drag is a minimum. Hence each of the lineal distributions of lift is 
elliptic and each contains the whole lift of the wing. Therefore the drags of these 
distributions are inversely proportional to the squares* of their lengths. Also, the 
Appendix shows that all sections of the volume distribution of minimum drag are 
bi-convex parabolic. Hence each of the lineal distributions of sources is that of a 
Sears-Haack body equal in volume to the wing volume. Therefore the drags of the 
distributions are inversely proportional to the fourth powers of their lengths. 


After this brief discussion of the application of supersonic area rule to this 
case, we turn to Figs. 3 and 4, where the variation of K, and K,, with Mach number, 
slew angle and semi-span-to-length ratio is displayed. As we should expect from the 


*The exponent follows from a dimensional argument: 
2 D L 

D (*) =[I?] and so k=[/-*) 
q q 
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FiGuRE 2. Oblique cuts. 


discussion, K, and K,, vary in a similar fashion but K, varies the more rapidly, 
since the drags of the lineal distributions depend on a higher power of the length. 
As the Mach number tends to one, K, and K,, both tend to one for all span-to-length 
ratios and slew angles, since the factors are based on optima for this limiting case. 

Referring to Fig. 2(c), we see that, for unslewed wings, an increase in span/ 
length ratio at a fixed Mach number increases the length of the lineal distributions 
(compare A,B, with A,B,). Similarly an increase in Mach number at a fixed 
span-to-length ratio also increases the length of the lineal distributions (compare 


A,B, with A,B,). Thus in each case K, and K, are reduced, as can be seen in 
Figs. 3(a) and 4(a). 
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Ss 
0-20 
0°25 
0-35 
0°45 
12 16 18.20 22 24 16 18, 20 22 24 
(a) SLEW ANGLE 0 various s/£ (b) SLEW ANGLE 15° various s/f 
22r 
03 bs° | 


2 14 1 20 22 24 14 16 18 20 22 24 
M M 
\C) SLEW ANGLE VARIOUS S/¢ (d)S/ = 0-25,VARIOUS ANGLES 
OF SLEW 
FiGurE 3. Variation of K, with Mach number. 
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FiGuRE 4. Variation of K, with Mach number. 


For the slewed wings the effect of Mach number is more complicated. Clearly, 
if s/1 is near enough to 0-5 the wing is nearly circular, slew has little effect and K, 
and K,, still decrease with M. For a fixed angle of slew, v, not exceeding 45°, it is 
clear that s// cannot be less than 


(s/D min tan 


~ 


At values of s// near this lower bound, the values of K, and K,, rise as M increases 
above one, the variation of K, being again more extreme. Referring to Fig. 2(d), 
we see how, as M increases above 1, the longest lineal distribution increases in 
length from A,B, to A,B, but the shortest distribution decreases in length from 
A,B, to A;B,. Owing to the inverse dependence on the length, the effect of the 
shortening on K, and K, preponderates and they increase. For high enough Mach 
numbers the shortest distribution again begins to increase in length with M. Before 
this Mach number is reached, K, and K, will have started to fall. This fall can be 
seen to start at Mach numbers below 2 for Y=15°, s/1=0-2; Y=20°, s/1=0°25; 
v¥=30°, s/1=0°35 in Figs. 3 and 4. 
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Between M=1 and M =2, a wide range of values of K, and Ky is found. Ata 
fixed Mach number and slew angle the drag factors increase rapidly as s/J is 
reduced towards ($/J)min. In the limit the wing becomes a slewed lifting line, the 
diagonal of the rectangular 2s x / box which contains it; p, the ratio of its area to 
that of the box, becomes zero. Thus decreasing values of p are associated with 
increasing values of K, and K,. This phenomenon depends only on the linearised 
theory of inviscid attached flow and is distinct from the penalty in drag factors 
incurred in using a plan form of smaller p in a real flow. This penalty arises because 
the volume and load distributions needed to obtain low values of the drag factors 
on plan forms of low p tend to be incompatible with the behaviour of a fluid of 
non-vanishing viscosity. 


3.2. Lirr/DraG RATIO 


The lift/drag ratio obtained with the assumptions set out in Section 2 is now 
considered. For consistency with the calculations for delta-like plan forms in 
Ref. 1, we choose a Mach number of two and a volume parameter, r=V/S*/?, of 
0:04. The assumptions about skin friction are similar to those of Ref. 1. Asa 
first approximation it is supposed that structure weight depends primarily on the 
structural aspect ratio, i.e. the largest aspect ratio of a wing of the family obtained 
by varying the slew or sweepback of the plan form. It is, therefore, realistic to 
consider values of lift/drag ratio at cruise for wings of the same structural aspect 
ratio or, as is more convenient here, the same axis ratio b/a. 


Figure 5 shows the variation of (L/D).-, obtained as described in Section 2, 
with the geometrical slenderness ratio s// for various values of the axis ratio b/a. 
The curves terminate on the left at the value of s// corresponding to the unslewed 
wing (s//=a/2b). Those for b/a=>4 show clear maxima at values of s/1 between 
0:16 and 0-18, the maximum being sharper and larger for the larger values of b/a. 
Increasing b/a at constant s// means reducing p; for values of s/1 below about 0°26 
this results in an increase in L/D but, for values of s// above this, L/D falls as b/a 
is increased. This shows that, for higher values of s/l, the increases in the drag 
factors discussed in Section 3.1 outweigh the reduction in p. This is of little interest 
in the present case as the lift/drag ratios are low at these values of s/l. 


The more significant part of Fig. 5 is enlarged in Fig. 6 and level curves of 
slew angle v and plan form parameter p are added. It is now seen that the wings 
with the largest L/D for given axis ratio (=> 4) have slew angles between 15° and 
20°. The unslewed wings have p==/4, so the left-hand boundary of the curves 
for fixed b/a is the level curve for both Y=0 and p==/4. As mentioned already, 
increasing b/a at fixed s// corresponds to a reduction in p. The higher values of 
L/D can now be seen to correspond to values of p which are very low by compari- 
son with delta-like wings (p=4 for a delta wing). Since, in general terms, the 
smaller the value of p the more difficult it is to achieve low drag factors on a wing 
in a real flow, it is of interest to look at the variation of L/D at constant p. It is 
found that the optimum values of s// for constant p are a little lower than those for 
constant b/a (Fig. 6): about 0:13 to 0-14 instead of 0:16 to 0°18. 
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M=2, r=0-04). 
6 
2 + 4 


The curve for p=0°5 from Fig. 6 is superimposed, in Fig. 7, on Fig. 15 of 
Ref. 1. This shows cruising lift/drag ratio for slender delta-like wings with 
7 =0:04 at M=2, comparable skin friction, realistic values of K, and various values 
of K, and K,. The comparison shows that the aerodynamic efficiencies of slewed 
ellipses and delta-like plan forms are comparable at M =2, for similar values of p. 
It is known that a real flow resembling that postulated in the theoretical treatment 
given earlier will only be possible if p is large enough. In the absence of a thorough 
consideration of the real flow past a slewed elliptic wing, it is supposed that the 
lower bound of these values of p will be similar to that for a delta-like plan form 
and so the comparison is made at the same value of p. The lift/drag ratio for the 
slewed ellipses of Fig. 7 is greatest near s//=0-15. Fig. 6 shows that this corres- 
ponds to a slew angle a little below 15°, so that K,, is about 1-1 (Fig. 4) and K, is 
about 2 (Fig. 3). Since K, for the delta wings“ is about 0:8 and the skin friction 
drags are comparable, the difference between the best delta wings in Fig. 7 and the 
slewed ellipses is due primarily to the different volume-dependent wave-drags. 
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Figure 6. Lift/drag ratio for M=2, r=0°04 with 
(i) various axis ratios b/a (——————) 
(ii) various slew angles ¥ (— - — - —) 
(iii) various plan form parameters p ( -—- 
T — T 
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K =| 
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/ / Vi 
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/ w 
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vi J 
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Figure 7. Lift/drag ratios: M=2, r=0-04, p=0°5, slewed ellipses and deltas. 
—— p=0'5 from Fig. 6. 
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FiGure 8. Effect of thickness and Mach number on lift/drag ratio for b/a=10. 


At higher Mach numbers the maximum lift/drag ratio falls. Fig. 8 shows the 
variation of L/D for b/a=10, 7=0:04 for M=2 and M=5S, plotted here against 
the aerodynamic slenderness parameter, 8s/1. The peak is lower and less marked 
at M=5. It occurs at a higher value of £s/1, though at a lower value of s/J. The 
major axis of the ellipse becomes sonic at 8s/1=0-5, so a large part of the leading 
edge of the optimum ellipse for M=S is supersonic. 


Figure 8 also shows the effect of an increase in 7 from 0-04 to 0:06 at M =2. 
The maximum lift/drag ratio is reduced and occurs for a slightly more slender 
wing. It should be pointed out that, from equation (7), a wing of 6,000 ft.? plan 
area has a maximum thickness of only 6:2 ft. when 7 =0-04 and 9-3 ft. when 7 =0-06. 
One would expect the smaller value of + to give enough volume for a long-range 


transport aircraft, so the volume distribution which leads to the lowest drag for 
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given volume is not efficient for passenger-carriage. The method of the present 
paper does not allow us to estimate the effect of varying the wing thickness 
distribution or of adding a separate fuselage. 


4. Conclusions 


(a) At a Mach number of two, slewed elliptic wings with optimum distributions 
of thickness and lift have maximum lift/drag ratios for a given structural aspect 
ratio at angles of slew between 15° and 20° and at semi-span/length ratios of 0°16 
to 0:18. The slenderness ratio is thus a little less than has been found for delta-like 
plan forms. 


(b) Also at a Mach number of two, the predicted level of lift/drag ratio using 
similar assumptions is very much the same for slewed elliptic wings and for slender 
delta-like wings. 
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Appendix 


THE SLEWED ELLiptTic WING OF LEAST DRAG FOR GIVEN VOLUME 


According to the supersonic area rule (see Ref. 5, for example), the drag due to 
volume of a thin wing is the average of the drags of the oblique equivalent bodies of 
revolution. Suppose the wing is cut by a family of parallel Mach planes which cut the 
wing plane in lines making an angle cot-'(8cos@) (-zr<=6< 7) with the flight 
direction, and that the sections so formed are projected on to planes normal to the flight 
direction. Then the streamwise distribution of these projected areas (called the oblique 
area distribution) is that of the oblique equivalent body of revolution corresponding to 
the angle 6. If the Mach plane is replaced by a plane normal to the wing plane passing 
through the intersection of the wing plane and the Mach plane, the difference between 
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the projected areas of the resulting cross sections is of the order of quantities neglected 
in thin-wing theory. Since this change simplifies the geometrical argument, it will be 
made in what follows. 


Each oblique equivalent body has the same volume as the wing. For a wing of 
elliptic plan form with a continuously curved surface, the oblique area distributions have 
continuous first derivatives which vanish at the end-points; therefore the drags of the 
equivalent bodies are least for given volume when they have Sears-Haack area distribu- 
tions. Thus the slewed elliptic wing would have least drag for given volume if all its 
oblique area distributions could be of the Sears-Haack type. 


The Sears-Haack distribution is the product of an elliptic distribution and a 
parabolic distribution. The lengths of the chords cut by parallel planes on an elliptic 
plan form are distributed elliptically. Thus, if all sections of the wing were parabolic, 
each oblique area distribution would be the product of an elliptic and a parabolic 
distribution and the drag of the wing would be a minimum for given volume. If the 
upper surface of the wing is a paraboloid (i.e. a quadric surface which touches the plane 
at infinity) whose axial direction is normal to the wing plane, all the sections are 
parabolic and the drag is a minimum. 


In this Appendix the drag of the wing so defined is calculated. In the special case 
of zero slew, the result agrees with that of Jones“? when an obvious correction is made 


to equation (71) of Ref. 2. We use the notation of Ref. 2, which is made clear in Fig. 1. 
The equation of the plan form is 


F (x, y)=x?—2m’ xy+(m? +a’ y?—a’2=90, . 
where b’=(a* cos? ¥+b? sin? y)!/? 
m’=(b*—a*) cos v sin 
a’=ab/b’. 
a is the semi-minor and b the semi-major axis of the ellipse and wv is the angle of slew 
between the major axis and the stream direction, Ox. If every plane normal to the wing 
plane, z=0, cuts the upper surface of the wing in a parabolic arc and the plane z=0 


cuts it in F (x, y)=0, the equation to the upper surface is z—AF (x, y)=0. A is specified 
in terms of the maximum thickness, f, which occurs at x=y=0, as follows :— 


t=2z=2A F (0, 0)= —2Aa”. 


The upper surface is therefore 
Z= | 2m’xy— x? + a/b") y? +a]. . 


The plane +my, (15) 
where m= cos 6, cuts the plan form where 


—b? (m—m’) x, + (d?—x,")] 


y= 
d? 


(16) 
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[hese points coincide when x,= +d. The two corresponding planes touch the plan 
form and so the length of the oblique area distribution is 2d. The length of the chord 
cut by one of the family of planes is, by equation (16), 


2d b’ 


(d?—x,”) sec (tan-! m). . (18) 


Each of the planes cuts the upper surface in a parabola, so the maximum thickness of 
the section is at the mid-point of the chord. The mid-points of parallel chords of an 
ellipse lie on a diameter and the section by the normal plane through a diameter is a 
parabola through the point of maximum thickness of the wing. Hence the thickness 
distribution of this section through a diameter is 


and this must be the maximum thickness of the section cut by the plane (15). 


The area of a bi-convex parabolic section of length L and thickness T is 2LT/3 
and so the projection on to a plane x=constant of the area cut by (15) on the wing is, 
using equations (18) and (19), 


4 abt x,*\ 
3d 


This is the oblique area distribution corresponding to the angle 4 (note: d=d(@)). It is 
a Sears-Haack distribution corresponding to a volume 


a maximum cross-sectional area 4abt/3d, a length 2d and a drag D (4), where 


D (@) 
=2na*b?1?/d*. . (21) 
q 
The drag of the wing is therefore 
D 
= dé 
q 2n q 
2s 
=a" (22) 


| [a’? +b’? (8 cos6—m’ffP 
To evaluate this integral, we consider the related integral 


dé 


J (a, b’, m’, * 
J (8 cos @—m’)* 


We also consider the remark that the integrand of J is the real part of 


d a +ib’ (8B cos 


4 
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Now, except when B/A is real and greater than one, 


dé 
A+B cos 6 (A?—B?) 


2x | 
Therefore \a’ J [a —im’ + 
ay 2a’ dé 
0a J) (8 cos 
D 1 oJ 
so that =a*b*t? (- 
q -ad ca 


, 


1 


If the ellipse is not slewed, Y=0, b’=a, m’=0, a’ =b and so 


_ D B*+2 (b?/a*) 
mabq 


(e.g. p. 24 of Ref. 6). 


(24) 


For 8=1, and a and b interchanged, (24) is the function plotted in Fig. 8 of Ref. 2. 
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The Elastic Stress Concentration Factors 
in Shouldered Shafts 


PART II: Shafts Subjected to Bending 


I. M. ALLISON 


(Department of Civil and Municipal Engineering, University College, London) 


4. Introduction to Part I 


This paper, which describes a photoelastic determination of the elastic stress 
concentration factors associated with shouldered shafts subjected to pure bending, 
forms the second part of an investigation initiated by the Structures Committee of 
the Royal Aeronautical Society. A general description of the manufacture of the 
models, the photoelastic technique employed to evaluate the s.c.f. and the results of 
the tests on shouldered shafts subjected to torsion have been given in Part I (Ref. 1). 


NOTATION (see Fig. 8) 

D large diameter of shaft 
d small diameter of shaft 
r shoulder radius 
shoulder depth (thus D = d+ 2h) 

r,@,2 cylindrical polar co-ordinates 
z axis of shaft 

P,Q,R_ principal stresses 

p.q.? directions of principal stresses 
n direction of the normal to the shaft surface 
ss normal stress tangential to surface in a direction s. 


The stress concentration factors are based on the nominal stress in the smaller 
diameter of the shouldered shaft. The factor K used throughout this paper is 
therefore defined as 


peak stress in the component 


maximum stress in a shaft of the 
smaller diameter under the same loading conditions 


K 


Received February 1961. Part I was published in May 1961. 
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TABLE | 
0-9 0-8 0-4 0-2 0-05 


1-0 | 8 in. long 
2-0 models 
3-0 NP 

4-0 NP 
5 
6 


67 | 10 1n. long 
models 
15-0 | 
95-0 


NP photoelastic test not practicable. 


The dimensions of the sixteen photoelastic models tested are specified by the 
parameters in Table I. 


Pure bending was applied to each of the standard photoelastic models, in the 
loading frame used for the torsion tests, by means of a suitable arrangement of links 
and wires, and the loaded Araldite model was subjected to an annealing cycle. 
Slices were cut from the frozen model and examined in a standard polariscope. 


Bending Tests 


-L.) STRESS PATTERN IN PURE BENDING 


The axial plane containing the peak stress is immediately apparent in the 
strained model, the peak stress being a principal stress (ss) lying in the axial plane 
and located on the surface of the shaft (see Fig. 8). The principal stresses at the 
opposite ends of a diameter are equal in magnitude but opposite in sign; thus, taking 
the numerical mean of the extrapolated results obtained from each of the axial slice 
boundaries eliminates the effects of time edge stress and any departures from the 
pure bending condition which might be induced by the loading system. 


5.2. PrEaK STRESS MEASUREMENT 


The peak stress was obtained from a slice taken in an axial plane, as in the 
torsion test series’?. The preparation of the slice and measurement of the relative 
retardation at points along the line OC remained the same. except that the slice was 
viewed at normal incidence instead of obliquely at angles of incidence at +45 and 

45°. Measurements at points along the lines OC (shown in Fig. 8) were made at 
both boundaries of the slice and the numerical mean of the extrapolated values of 
the relative retardation at the surface was obtained. Since the tangential surface 
stresses (ss) on opposite boundaries are equal in magnitude but opposite in sign, the 
averaging process immediately eliminates the following five possible sources 
of error. 
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(i) Any tangential stress associated with initial hoop strain, induced 
by the thermal gradients which are present during casting or 
curing of the Araldite resin. 

(ii) Axial stress due to friction. 

(iii) The effects of self weight. 
(iv) “Edge stress,” a tangential compression produced by the 
absorption of water from the atmosphere. 

(v) Secondary effects due to any change in shape of the shoulder 
root radius associated with the relatively high applied strain 


(approximately 2-5 per cent maximum in the frozen model). 


5.3. “UNIFORM” STRESS MEASUREMENT 


The “uniform” stress was obtained from a slice cut in the same axial plane as 
that of the “peak stress” and located away from the disturbing effect of the shoulder. 
Observations of the relative retardation at normal incidence were made at points 
along a diameter in the smaller shaft. The results were plotted and extrapolated to 
obtain the tangential stress at the surface. For each test the linear relationship 
between the relative retardation and the distance from the shaft surface was checked 
in order to confirm that the model behaved elastically and that the applied load was 
sensibly a pure bending moment. This linear behaviour also confirmed that the 
chosen section was sufficiently far from the disturbing effect of the shoulder and that 
the level of initial stress in the model was negligible. The numerical mean of the 
theoretically equal but opposite tangential stresses at the surface was taken in order 
to eliminate the effects of time edge stress or any small axial load which might have 
been induced in the model. 


The measurement of the nominal stress in a slice cut from the model itself offers 
two further advantages in that it eliminates both the calibration of the photoelastic 
material and the need to compute the magnitude of the applied bending moment. 
In a photoelastic test the usual practice is to calibrate the photoelastic material by 
making measurements on a tension test piece manufactured from a plate which has 
been poured from the same mix as the model under test. However, even with resin 
of the same composition, identical optical properties cannot be guaranteed for 
sections which are as different as a thin plate and a cylindrical shaft. 


Ficure 8. Notation and slices for the bending tests. 
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ASYMPTOTE CONDITION 
| FOR h/D=0°4 | 


FiGURE 9. Stress concentration factor K against h/r, for constant values of h/D. 
Load: pure bending. 


Also, a dead weight loading system of the type used to stress the models is 
inevitably subject to friction, particularly at high temperatures (120°C); thus the 
accurate measurement of the actual moment applied to the model presents 
difficulties. In the present series of tests both these problems have been avoided by 
the method used to evaluate the nominal stress. 


5.4. RESULTS 


The stress concentration factor is defined as 


peak stress in the shouldered shaft 
maximum axial stress in a uniform shaft of the ° 
smaller diameter, subjected to the same moment 
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FiGurE 10. Stress concentration 
factor K against d/D. for con- 


stant values of r/d, Load: pure ra 
bending. 
4 

/ 

1 


As in the torsion tests on shouldered shafts the families of curves for constant values 
of h/ D (Fig. 9) and constant values of r/d (Fig. 10) have been selected from a series 
of plots of K against different dimension parameters. 


The general shape of the families of curves for the bending tests is similar in 
form to that obtained in the other shouldered shaft tests. Four points of similarity 
to the torsion test results are immediately apparent, as follows : — 


(i) The s.c.f. becomes very large as r/d tends to 0. 
(ii) For all values of r/d the s.c.f. tends to 1 as d/D tends to 1. 
(iii) For all values of h/D the s.c.f. tends to 1 as A/r tends to 1. 


(The effects (ii) and (iii) would be expected from simple reasoning). 


(iv) As the value of d/D is decreased from 1 the s.c.f. increases 
initially, reaches a maximum at d/D ~ 0-65, and then decreases. 
This behaviour is apparent for all values of r/d, but is more 
marked for the smaller values of this parameter. As in the 
torsion tests, it is reasonable to expect the s.c.f. to increase as 
d/D is decreased from 1, (i.e. as the shoulder depth A is 
increased) but intuitive reasoning would predict that a decrease 
in the ratio d/D by the addition of extra material to an already 
large shoulder would make negligible difference to the s.c.f. 
Clearly any simple argument such as that based on the two- 
dimensional behaviour in flat shouldered plates, or that 
postulating the existence of the plane strain condition in the 
shaft, will inevitably lead to a solution in which K increases and 
becomes asymptotic to a limit as d/D is decreased to zero. 
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Predictions of the qualitative behaviour of a stress raiser in a complicated 
system by the extension of arguments which only strictly apply to a simple case are 
widespread and often indicate the correct trend. Before discussing the reasons for 
the decrease in K it is necessary therefore to establish the significance of the 
experimental curves as drawn by the author, particularly when compared with the 
assumption that the results should become asymptotic as d/D tends to zero. 


5.5. THE SIGNIFICANCE OF THE EXPERIMENTAL RESULTS 


If it is assumed that the s.c.f. for each value of r/d increases monotonically and 
becomes asymptotic to the maximum value obtained in the photoelastic tests (e.g. 
for r/d=0-1, Fig. 10 gives the asymptote as K =1-:75), but shows no subsequent 
decrease as d/D decreases to zero, then the effect is to raise the level of the lower 
curves h/D=constant in Fig. 9. In particular, the curves h/ D=0°3 and h/ D=0-4 
which satisfy the asymptote condition are shown by chain-dotted lines in Fig. 9. 
The difference between each of the experimental results and the appropriate value 
of the asymptote curve has been expressed as a percentage of the asymptotic value. 
This difference is described in Table II as the percentage error on an assumed 
behaviour “assumes. A Comparison of these values with the corresponding 
difference ““foxperimentar’» Detween the experimental results and the author’s curves 
(shown as full lines in Fig. 9) is given in Table II. 


100 (Ky, — Ks) 100 (K,, K) 


where K,=s.c.f. obtained from photoelastic test. 
K =corresponding s.c.f. obtained from the curves of “best fit” drawn 
through the experimental results. 
K,s=corresponding s.c.f. obtained from the curves satisfying the 
asymptote condition. 
Two points are immediately apparent from an inspection of Table II: 

(i) The error in the experimental curves is random and never 
exceeds 7 per cent. 

(ii) The error between the results of the tests and the asymptote 
curves is always of the same sign and, for all but the smallest 
value of h/r, where the difference is in any case negligible. is 
greater than the maximum estimated error in the experimental 
points. 


TABLE II 
A COMPARISON OF THE PERCENTAGE DIFFERENCES FOR THE EXPERIMENTAL CURVES AND THE CURVES 
SATISFYING THE “‘ASYMPTOTE”’ CONDITION 


h/D 0:3 0-3 0-4 0-4 0-4 0-4 0:4 0-475 

h/r 7:5 15 1-0 3-0 4-0 10 16 95 
Esssumed 8 -1°5 -9 -18 -24 38 
Cexper mental 0 6 4 +2 9 0 


-6°8 -5 0 
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Although these facts were established from the results of the tests included in 
the original programme, a relatively small number of experimental results was 
available in the region where the decrease in s.c.f. became marked. Also, the 
parameters specifying the dimensions of the “insert” model (h/D=0-475, h/r=95) 
were so different from those of the bulk of the models tested that the result of the 
insert test was of little assistance in determining the shape of the curves in the region 
0:4>d/D> 0-2. However, this result did provide confirmation of the fact that 
K decreased for small values of d/D. 


The lower limit of the parameter d/D to be included in the original programme 
was decided from practical considerations of the maximum increase in shaft 
diameter which was likely to occur in practice. In order to confirm that marked 
differences between the asymptote values and experimental results could occur, it 
was decided to undertake two further tests on models having dimensions outside 
the range normally encountered. The parameters specifying the dimensions of 
these models (h/D=0-4, h/r=10 and h/D=0-4, h/r=16) were chosen to indicate 
the maximum divergence from a single, reasonably small, diameter casting. The 
final large diameters of the two models were 1} in. and 2 in. respectively. The 
results of these tests are plotted in Fig. 9 and the relevant differences are given in 
Table II. The experimental results provide adequate confirmation of the original 
tests. The shape of the asymptote curves shown in Fig. 9 cannot be reconciled with 
the experimental points within the limits of the maximum possible estimated error; 
thus the results of the present tests must predict a decrease in the stress concentration 
factor for larger values of d/D. 


It is interesting to note that the results of the two additional tests lie within 
5 per cent of the values estimated by extrapolating the curves drawn through the 
original experimental results. The close agreement between the results of the 
additional tests and this extensive extrapolation confirms the increased accuracy 
inherent in families of curves which have been drawn to satisfy five independent 
cross-plotting conditions. 


5.6. POSSIBLE REASONS FOR THE DECREASE IN K 


The boundary of the shouldered shaft is, in general, composed of four distinct 
regions, as follows : — 
the cylindrical portion of small diameter, 
the blend radius at the root of the shoulder. 
the plane face of the shoulder (which decreases and finally disappears 
as the blend radius is increased), 
the cylindrical portion of large diameter. 


The author knows of no theoretical solution for the shouldered shaft when 
subjected to torsion or bending. The absence of an exact solution is probably due 
to the difficulty of satisfying this complicated boundary. However, there is no 
theoretical reason why the s.c.f. should not increase initially and then decrease as 
d/D is decreased, and it can be seen from Fig. 9 that whether the decrease occurs 
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Predictions of the qualitative behaviour of a stress raiser in a complicated 
system by the extension of arguments which only strictly apply to a simple case are 
widespread and often indicate the correct trend. Before discussing the reasons for 
the decrease in K it is necessary therefore to establish the significance of the 
experimental curves as drawn by the author, particularly when compared with the 
assumption that the results should become asymptotic as d/D tends to zero. 


5.5. THE SIGNIFICANCE OF THE EXPERIMENTAL RESULTS 


If it is assumed that the s.c.f. for each value of r/d increases monotonically and 
becomes asymptotic to the maximum value obtained in the photoelastic tests (e.x. 
for r/d=0-1, Fig. 10 gives the asymptote as K = 1-75), but shows no subsequent 
decrease as d/D decreases to zero, then the effect is to raise the level of the lower 
curves //D=constant in Fig. 9. In particular, the curves h/D=0°3 and h/ D=0°4 
which satisfy the asymptote condition are shown by chain-dotted lines in Fig. 9. 
The difference between each of the experimental results and the appropriate value 
of the asymptote curve has been expressed as a percentage of the asymptotic value. 
This difference is described in Table II as the percentage error on an assumed 
behaviour “assumes A Comparison of these values with the corresponding 
difference “experimental, between the experimental results and the author’s curves 
(shown as full lines in Fig. 9) is given in Table IIL. 


100 (Ky — 100 (K,, — K) 


where Ky =s.c.f. obtained from photoelastic test. 
K =corresponding s.c.f. obtained from the curves of “best fit” drawn 
through the experimental results. 
K,s=corresponding s.c.f. obtained from the curves satisfying the 
asymptote condition. 
Two points are immediately apparent from an inspection of Table II: 

(i) The error in the experimental curves is random and never 
exceeds 7 per cent. 

(ii) The error between the results of the tests and the asymptote 
curves is always of the same sign and, for all but the smallest 
value of h/r, where the difference is in any case negligible, is 
greater than the maximum estimated error in the experimental 
points. 


TABLE II 
A COMPARISON OF THE PERCENTAGE DIFFERENCES FOR THE EXPERIMENTAL CURVES AND THE CURVES 
SATISFYING THE “‘ASYMPTOTE” CONDITION 


h/D 0-3 0-3 0:4 0-4 0:4 0:4 0:4 0-475 
h/r 7-5 15 a) 3-0 4-0 10 16 95 
“assumed 8 -13 ~9 -18 —24 -23 38 
= x} mental 0-6 4 +2-9 +0-8 0 
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Although these facts were established from the results of the tests included in 
the original programme, a relatively small number of experimental results was 
available in the region where the decrease in s.c.f. became marked. Also. the 
parameters specifying the dimensions of the “insert” model (h/D=0-475, h/r=95) 
were so different from those of the bulk of the models tested that the result of the 
insert test was of little assistance in determining the shape of the curves in the region 
0:4>d/D> 0-2. However, this result did provide confirmation of the fact that 
K decreased for small values of d/D. 


The lower limit of the parameter d/D to be included in the original programme 
was decided from practical considerations of the maximum increase in shaft 
diameter which was likely to occur in practice. In order to confirm that marked 
differences between the asymptote values and experimental results could occur, it 
was decided to undertake two further tests on models having dimensions outside 
the range normally encountered. The parameters specifying the dimensions of 
these models (h/D=0-4, h/r=10 and h/D=0-4, h/r=16) were chosen to indicate 
the maximum divergence from a single, reasonably small, diameter casting. The 
final large diameters of the two models were 1} in. and 2 in. respectively. The 
results of these tests are plotted in Fig. 9 and the relevant differences are given in 
Table II. The experimental results provide adequate confirmation of the original 
tests. The shape of the asymptote curves shown in Fig. 9 cannot be reconciled with 
the experimental points within the limits of the maximum possible estimated error; 
thus the results of the present tests must predict a decrease in the stress concentration 
factor for larger values of d/D. 


It is interesting to note that the results of the two additional tests lie within 
5 per cent of the values estimated by extrapolating the curves drawn through the 
original experimental results. The close agreement between the results of the 
additional tests and this extensive extrapolation confirms the increased accuracy 
inherent in families of curves which have been drawn to satisfy five independent 
cross-plotting conditions. 


5.6. POSSIBLE REASONS FOR THE DECREASE IN K 


The boundary of the shouldered shaft is, in general, composed of four distinct 
regions, as follows : — 
the cylindrical portion of small diameter, 
the blend radius at the root of the shoulder. 
the plane face of the shoulder (which decreases and finally disappears 
as the blend radius is increased), 
the cylindrical portion of large diameter. 


The author knows of no theoretical solution for the shouldered shaft when 
subjected to torsion or bending. The absence of an exact solution is probably due 
to the difficulty of satisfying this complicated boundary. However, there is no 
theoretical reason why the s.c.f. should not increase initially and then decrease as 
d/D is decreased, and it can be seen from Fig. 9 that whether the decrease occurs 
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or not is simply dependent on the relative rates of increase of K with h/r, for 
constant values of h/D. 


It is unlikely that the s.c.f. is a simple function of d/D and r/d only, and a 
combination of two separate effects is likely to determine the value of K for any 
particular value of d/D in Fig. 10. These two effects can be described qualitatively 
as follows: - 

(i) The diffusion of the axial stress around the blend radius. 
(ii) The rate at which the general level of stress decreases in passing 
from the small diameter shaft to the larger in the axial direction. 


For the larger values of d/D, i.e. small shoulder depths, the first effect will be 
predominant, since the general level of stress in the small and large diameter shafts 
does not change excessively, and the axial stress will diffuse almost completely 
around the blend radius. The axial stress on the surface of the large diameter shaft 
increases from zero, at the shoulder, to its maximum value in a fairly short length 
of the shaft and the value of K will be almost entirely a function of the blend radius 
alone. As d/D is decreased the second effect will become increasingly important, 
especially as the decrease in stress level at a particular radius, for bending or torsion, 
is proportional to (d/D)'. It is possible that in the presence of large shoulders the 
second factor inhibits the diffusion of stress around the blend radius to such an 
extent that a decrease in the s.c.f. is produced. Under these conditions the s.c.f., for 
a particular value of r/d, would increase initially as d/D is decreased, but would 
ultimately reach a limit and then decrease at a continuously increasing rate as the 
“inhibiting” factor becomes dominant. 


5.7. COMPARISON WITH THE RESULTS OF PREVIOUS INVESTIGATIONS 


Three other investigators have obtained experimental values for the s.c.f. in a 
shouldered shaft subjected to bending. Their results are compared with those of 
the present investigation in Fig. 10. 


Peterson and Wahl? made measurements with a mechanical strain gauge on a 
large steel shaft having a D/d ratio of 1:5 and obtained results for three values of 
r/d (namely 0:17, 0:28 and 0:5). It can be seen that these results lie within 8 per 
cent of those of the present investigation. Thum and Bautz employed a similar 
method and obtained the curve, shown chain-dotted in Fig. 11, for a steel shaft also 
having a D/d ratio of 1:5. This curve is in reasonable agreement with the results 
of the present investigation over the practical range of shoulder blend radii. 
However, the results of Thum and Bautz do not predict a rapid increase in the s.c.f. 
for values of r/d less than 0:06. This divergence in the experimental results for low 
values of r/d is to be expected, since the use of a mechanical extensometer is likely 
to predict low values of peak stress for very small blend radii. Also the accuracy 
of measurement for both experimental methods decreases rapidly at the lower end 
of the r/d scale. 


Using the photoelastic stress freezing method, Leven and Hartman have 
obtained results with three-dimensional models for D/d ratios of 1:25 and 2-0. All 
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Ficure 11. Comparison of the results of the present investigation with those of previous tests. 


the experimental points, with the exception of two (specified by the parameters 
D/d=1:25, r/d=0-07 and D/d=1-25, r/d=0-063), lie within 5 per cent of the 
curves obtained in the present series of tests. An extensive extrapolation of these 
curves indicates that the results of the two photoelastic investigations are not incon- 
sistent even for r/d values as low as 0-020, although the actual value of the s.c.f. in 
this region must be subject to a wide error band in view of the inaccuracy incurred 
by extrapolating curves of this type. 


The comparison of the results with those of previous investigators, and an 
estimation of the probable errors in the photoelastic measurements, suggests that 
the maximum error in the curves in the range 0:05 << r/d<1,04<d/D<1 is 
not greater than +6 per cent. 
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A Comparison of the Characteristic Equations 
in the Theory of Circular Cylindrical Shells 


D. S. HOUGHTON, M.Sc.(Eng.) and D. J. JOHNS, M.Sc.(Eng.) 
(The College of Aeronautics, Cranfield) 


SUMMARY: Characteristic equations are derived for thin circular shells, based 

on various approximations to the linear elastic theory of small deformations. 

By representing the deformation in a Fourier series in the circumferential 

direction, the roots of these equations are computed for a range of the 
significant parameters and compared. 


1. Introduction 


In recent years, a number of authors have presented approximate solutions to 
the small deflection theory of circular cylindrical shells with particular emphasis on 
problems concerned with local loadings. In a paper by Jaeger and Chilver*’, 
various solutions have been examined and the corresponding characteristic equations 
presented. In addition, numerical solutions were presented to the equation derived 
by Biezeno and Grammel?. 

The purpose of this paper is to examine the characteristic equations obtained 
from various shell theories. It is not the intention to compare critically the basic 
analyses and assumptions of these various theories but, rather, to compare the form 
of, and the solution to, the resulting characteristic equations. 

It is believed that refinements of Love’s approximation, such as are inherent in 
many of the theories considered, are, in general, meaningless because of the neglect 
in all the theories of certain other effects. However, in view of Novozhilov’s“ own 
comments (see Section 3) with regard to the Love approximation, from which 
Timoshenko’s™ equations follow, these two sets of equations were deemed suitable 
to indicate the manner in which the characteristic equations are obtained from the 
three equilibrium equations. 


NOTATION 
a_ radius of cylinder 
A,B,C arbitrary constants in expressions for u, v, and w 
p number of circumferential waves 
Received February 1961. 
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t thickness of shell 
u,v,w displacement of a point in the median plane of shell 
x,y axial and circumferential co-ordinates 
a= ft? /(12a) 
(=x/a 
X exponential index in expressions for u,v and w 
» Poisson’s ratio 


o=y/a 


2. Theory 


2.1. GENERAI 


Love’? has presented general equations defining the linear elastic deformation 
of an element of a shell having arbitrary curvature, which have been developed by 
Timoshenko for the circular cylindrical shell. Similar equations have also been 
obtained by Novozhilov, who has discussed the assumptions made by Love and 
other authors, and who states that the development given by Love is not entirely 
free from inadequacies. However, examination of the equations shows that the 
differences in the two solutions are confined to small-order terms. 


In the following analyses, it is assumed that the inner and outer surfaces of the 
shell are free from loads. Hence, the equations developed by Timoshenko are 
given as 


Ox 2 oxdy aox 
l+v O7u C74 C*y C’w Cw 
2 dxdy ox? ay’ Ox*dy ady 
(1b) 
—— +2a |(2-v) + 35] + = + 2a’°V‘w=0. (le) 
aox acy Ox?dy ov a 


For the same cylindrical shell the Novozhilov equations have been derived in 
Ref. 6, from which it may be seen that they are identical to equations (1) except for 
equation (1b) which becomes 
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2 oexdy 2 ox? dy?’ ox* dy 
Ow 
+ aa + 33] =0. (1B) 
Ox?dy dy 2dy 
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It is obvious, from inspection of the terms in v in equations (15) and (1B), that 
for small values of 2 the corresponding terms in the two equations are negligible. 
This approximation has been made in Ref. 7. In fact, according to Vlasov‘, all 
of the square-bracketed terms in equations (1) are small quantities of the same order 
as those which have been disregarded in the derivation of these equations. This 
comment is not justifiable by inspection but, assuming it to be correct, the following 
equations are obtained from both equations (1) and Ref. 6. 


14+) 


aw 

2 dy? 2 dxdy~ adx 

l+v 1—vo7v ‘ Ow 0 (2b) 

2 oxdy 2 dy? ady ° 

vou ov w 

+—+2aaViw=0. . 2c 

aox ady @ (2c) 


It is well known that an eighth-order partial differential equation may be 
derived in any one of the displacements by elimination of the other two displace- 
ments. However, by a suitable representation of the displacements u, v and w in the 
form of a double Fourier series, the determinantal equation formed from the three 
equations in u, v, and w can be expanded to give the characteristic equation. 


In the present analysis the modes of deformation to be considered are 
= 
u= & A,e sin po 
) 


p=0 


w= C,e™ sin pp 
p=0 
These forms are particularly useful for a complete cylindrical shell carrying 
arbitrary loads on the circumferential boundaries. 


The condition for non-zero values of A,, B, and C, is that the determinant of 
the coefficients of A,, B, and C, should vanish on substitution of equations (3) into 
equations (1). 


2.2. TIMOSHENKO’S EQUATION 


Substitution of equations (3) into equations (1), and expansion of the 
determinantal equation so formed, gives the following characteristic equation if 2 is 
considerably less than unity. 


(3+) 


—2p? 307+ | pt (4) 
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If the 2-term in v in equation (1b) is assumed te be negligible, the corresponding 
characteristic equation becomes 


AS — + 6p* — 2p’ |20*- 507 + v) A? 4 

+ p® (p?—2)=0. (5) 
This equation was originally obtained by Bijlaard“ who made the foregoing 
assumption, and it is presented with equation (4) in Table I. 
2.3. NOVOZHILOV’S EQUATION 


In a similar manner to that in Section 2.2, Ref. 6 gives the following 
characteristic equation : — 


+ p* (p?-—1)?=0 (6) 


If the z-term in v in equation (1B) is assumed to be negligible, the corresponding 


result is 
y2 


+ 6p? (p? - ~ 4p* (p? —2) A? + p* (p? —2) 30. 


(7) 


Equations (6) and (7) are also given in Table I. 


2.4. DONNELL’S EQUATION 


The characteristic equation corresponding to equations (2) based on the 
simplifying assumption of Vlasov is 


2 


Equation (8) has also been derived from Donnell’s® equation in Ref. | and it is 
seen from Table I that, by comparison with the other solutions, only the higher 
order terms in p have been retained. Evidently the accuracy of equation (8) 
improves as p increases. 

A further approximation suggested by Biezeno and Grammel® may be 
made when 


1 —v? 
6p‘ > 
in which case equation (8) reduces to 
+ . . . 
or . . (10) 


and the solution is independent of 2. 
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Biezeno and 


Grammel'? 


Fliigge™ 


Vlasov 


Morley 


limoshenko" 


(equation (4)) 
Bijlaard'’ 

(equation (5)) 
Novozhilov™ 

(equation (6)) 


Equation (7) 


Donnell'® 


fequation (8)) 


Naghdi and 


Berry‘! 


Kennard'!*) 
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TABLE I 
COMPARISON OF GENERAL CHARACTERISTIC EQUATIONS 


| 
AS—2 (2p?—v) AS 4 | + \4— 


— 2p? (4—v) p?+(2— A? + p* (p?—-1P=0. 


{—v2 
AS—2 (2p?—v) + 


Ea — A2 + p* (p?—1P=0. 


4 6p? 


\ 


— 2p? ~ ) +1] A? + (p*—1=0 
T6p* (p*—1)4 | 3p? +1] \24 


+ p* (p?—1)?=0. 


E 
2 (2p? — 1) A* + 


x 


v) 
) A? + p* (p?—1)?=0. 


A8—2 (2p?) AS + 


—2p? | 2p'- | ) \? + p* —2p?)=0. 


—v2 
— 4p? + p* (p?—- 1 =0. 
1» 
= 2 (ap + + 6p* (p?—1)—2 (1 —v?) p? | A*— 4p? (pt —2p?) AF + 


+ p* —2p?)=0. 


+ Af- =p" (2p*) A2 4 p* 


—2p? | 2p*—(3+v) p?+ A? + p4 (p?— =0. 
| @ 2(1—v) 


(8 —8v — (2+ v)(2—3v 
) A* + p' 1)?=0. 


4(1—v) 
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TABLE I 
ROOTS OF CHARACTERISTIC EQUATIONS FOR p=0 


AXISYMMETRIC DEFORMATION MODES 


25x10 104 =2:5xX10 
x a 
Biezeno and 498685 9-99273 49-9985 
Grammel 501684 + 1000773: 50-0015: 
498498 999250 49-9985 
501498; 1000750; + 50-0015 
498548 9:99256 49-9985 
501548: 10-007 56i + 50-001 5: 
4:95025 9-97503 49-9950 

505025 1002503; 50-0050; 
Timoshenko!! 5:0009 1 10°00011 50-0000 
(equation (4)) 500091; 10-0001 1; + §0:0000i 
Biulaard 5.00000 10-00000 50-0000 
(equation (5)) 5-00000i + 10:00000; + 50-0000; 
Novozhilov 5-00182 10°00023 50°0000 
(equation (6)) + 500182 + 10°00023; + 50°0000; 
5-00000 10-00000 50-0000 

Equation (7) 
5-00000 10: 00000; + 50-0000; 
Donnell'* 500000 10:00000 50-0000 
(equation (8)) + $-00000/ 10-00000/ + $0°0000i 
Naghdi and 5-00000 10-00000 50-0000 
Berry ! + §-00000/ 10: 00000; + 50°0000% 
+ 500082 10:00010 50-0000 

Kennard 
+ 5-00082i + + §0-00007 


3. Presentation of Results 


In addition to the characteristic equations which are derived in Section 2, 
Table I also presents a number of other equations which are found in the 
literature (Refs. 10-13). Solutions of the equations in A for a range of the 
parameters (1 —v?)/2z and p have been derived and the results are given in Tables II 
and III. In these solutions, Poisson’s ratio has been taken as 0-3 and the values 
taken for (1 — v?)/ are those suggested by Hoff” and used by Jaeger and Chilver. 
When (1 —v?)2 < the shell is relatively thick, and the simplifying 
assumptions made in deriving the basic equations may invalidate this theory; the 
value (1 —v?)/2=2°5x 10" probably represents a reasonable lower limit to shell 
thickness for many aircraft problems. 
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4. Discussion of Results 


It can be clearly seen from Table I that the effect of differences between the 
various characteristic equations becomes small as p increases. 


The equation suggested by Fliigge”’ is exactly the same as that proposed by 
Biezeno and Grammel?, except for a small term in the coefficient of A* (see Table I). 
Obviously as the term (1 — v*)/z increases, and as p increases, the influence of this 
term must diminish. The effect of this can be seen in Tables II and III. 


The effect of the term 2¥A° in some of the characteristic equations given in 
Table I is less pronounced as p increases. In fact, for p=0 (Table ID), where the 
effect is greater, the roots of the various solutions do not appear to be too dissimilar, 
and their agreement improves as (1 — v*)/z increases. 


Morley’s solution” is based on an analysis which includes assumptions that 
are tantamount to assuming Poisson’s ratio to be unity, except in the term (1 —v?)/ 2, 
in the characteristic equation derived by Biezeno and Grammel. It is worth noting 
that, with equations such as those by Morley and Donnell, it is possible to determine 
the values of p without resort to the solution of a quartic equation, and it is seen 
from Table III that Morley’s solution is an improvement on that of Donnell, 
compared with the more accurate results of Fliigge“"’ and Biezeno and Grammel”. 

Hoff"* has suggested that Donnell equations should not be used when p << 4, 
and inspection of the results in Table III supports this view. However, for p=0 the 
results in Table II do not suggest that Donnell’s equation is significantly less 
accurate than the other more exact equations. 

Bijlaard™ has developed the shell equations put forward by Timoshenko™. The 
resulting characteristic equations ((4) and (5)) are very similar, and the agreement 
between the roots is excellent, particularly for large values of p and (1—v’*)/a. 
Similar observations can be made for equation (6), derived from Novozhilov’”, and 
equation (7) of this paper. 

The application of the characteristic equations to shell problems is not 
discussed in this paper, but reference can be made to Refs. 7, 15 and 16. 


5. Conclusions 

Characteristic equations are derived and presented for various approximate 
solutions to the theory of circular cylindrical shells. Solutions of these equations 
are given for a range of the parameters p and (i-—v’*)/2, and the roots of the 
characteristic equations are compared. 

In problems involving the use of the characteristic equations there is some 
advantage in using the special simplified equations such as those by Donnell and 
Morley, although Donnell’s equations should not be used for values of 0< p S3. 
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The Theory of Melting Ablation, with 
Vaporisation, Gas-Phase Chemical Reaction, 
Surface Pyrolysis, and Transient Effects 


D. B. SPALDING 


(*/mperial College of Science and Technology) 


SUMMARY: Standard techniques of mass-transfer theory are used for the 
prediction of ablation rates; the thermodynamic and chemical-kinetic material 
properties are introduced by way of enthalpy-composition diagrams; the inter- 
face condition is given as the root of a single non-linear equation involving 
material properties. The paper treats both steady ablation and unsteady 
ablation, the latter by means of a quasi-steady assumption which confines 
transient effects to the solid phase. 


Detailed comparisons are made with the methods of Sutton, Bethe and 
Adams, and Lees. The formulations of the paper are shown to be equivalent 
to, but more general than, those of the earlier authors; it is suggested that they 
are simpler as well. Some improvements over previous practice are 
recommended in connection with the calculation of the shear stress at 

the interface. 


1. Introduction 
1.1. THE PROBLEM CONSIDERED 


Missiles and satellites re-entering the earth’s atmosphere are subjected to very 
high rates of heating by the hot gases which are present in their boundary layers. 
Either by accident or design, these heat fluxes usually result in the removal of solid 
material from the missile nose; this material may either melt and flow along the 
surface, or it may vaporise or sublime. Sometimes the material vapour reacts 
chemically with gases in the boundary layer; and sometimes the change of phase 
itself involves a chemical reaction, commonly termed pyrolysis. The whole process 
of material removal is nowadays called ablation. 


It is obviously important, in the design of missiles, to be able to predict the 
amount and nature of the ablation process. How to make such predictions is the 
concern of the present paper. 


*This work was performed in the author's capacity as consultant to the Lockheed Missiles and 
Space Division, Palo Alto, California. 


Received January 1961. 
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1.2. OUTLINE OF PREVIOUS WoRK ON THE SUBJECT 


The ablation-prediction problem has been considered theoretically by several 
authors in recent years. Major contributions have been made by Sutton“'*: '”’, Bethe 
and Adams‘ and Lees“, as part of a large programme carried out in the United 
States; a paper by Adams“? gives a valuable survey of the achievements of this 
programme. 


However, ablation can also be regarded simply as a particular sort of mass- 
transfer process; it is thus amenable to analysis by means of techniques which have 
an earlier origin and which are in some respects more highly developed than those 
of the papers just referred to. This line of thought has been exploited by the present 
author®: '”, who has advocated the use of enthalpy-composition diagrams as aids 
to comprehension and calculation of ablation processes. Melting ablation was 
treated briefly in the first of these two references, but only for the case of very low 
liquid viscosity. 


1.3. PURPOSE OF THE PRESENT PAPER 


It is known from the work of Sutton“ and Bethe and Adams’, for example, 
that the viscosity of practically important materials is not “low” in the sense of the 
previous paragraph; these authors have shown how to calculate the effect on the 
ablation rate of finite and temperature-dependent liquid viscosities. One purpose of 
the present paper is therefore to incorporate this knowledge into the graphical 
procedures just referred to; the advantages of so doing are: (i) that calculation of 
ablation rates is thereby made easier, (i/) that the importance of the thermodynamic 
and chemical-kinetic properties of the materials is clarified, (iii) that the trends and 
interactions resulting from changed material properties or operating conditions can 
clearly be perceived, and (iv) that the comparison and choice of suitable nose-cone 
materials is facilitated. 


While the paper is primarily concerned with melting ablation, the opportunity 
will be taken to show the connections between the two modes of analysis for other 
types of ablation; in particular, terms such as “heat of ablation” and “heat- 
blocking” etc. will be related to concepts of conventional mass-transfer theory. It 
will also be shown that, in addition to the improvements in clarity of representation 
which results from the use of enthalpy-composition charts, detailed improvements 
can also be made to the theories of Bethe, Adams and others in respect of the 
solutions of the boundary-layer equations which they employ. 


Finally it is intended to show how the same techniques can be applied whether 
steady or unsteady ablation is in question. 


It might be argued that the treatments of ablation developed by Sutton, Bethe 
and Adams, and others, are quite satisfactory for the job they have had to do and 
that therefore there is no need for new approaches to the problem. If that were 
indeed the case, the present paper would have to be justified by pointing out a 
different purpose which it serves, namely that of a test of the general methods of 
mass-transfer analysis; for, if these methods are sufficiently general in concept, 
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convenient in notation, and flexible in manipulation, they should be able to 
accommodate the ablation problem without modification or extension. The present 
paper can be regarded as a study of whether this is indeed possible. 


1.4. OUTLINE OF THE PAPER AND ITS MAIN RESULTS 


Section 2 assembles the relevant standard results of mass-transfer theory and 
applies them to the ablation problem. It is shown that steady ablation rates can be 
calculated very simply in two important extreme cases; in the general case of 
melting ablation, however, uncertainty concerning the interface condition prevents 
immediate prediction of the ablation rate. The relation of the present formulation 
to those of other authors is discussed at some length in Section 2.3. 


In Section 3 we consider the hydrodynamics of the liquid layer and reduce the 
problem of determining the interface condition essentially to that of solving a single 
equation (equation (47)); this must usually be done numerically or graphically. 


Section 4 outlines procedures for predicting transient ablation phenomena. The 
analysis is a quasi-steady one; it reduces the problem to that of unsteady one- 
dimensional heat conduction with characteristic non-linear boundary conditions. 


NOTATION 


Note: Typical units are given, in which each quantity might be measured. 
(—) signifies “dimensionless.” The equations mentioned are those in which the 
symbol first appears. 


A dimensionless gas-phase drag coefficient (—) (equation (36)) 


B_ dimensionless driving force for mass transfer into the gas phase 
(—) (equation (1)) 


D dimensionless gas-phase mass-transfer conductance (—) (equation 
(45)) 


f mass fraction of nose-cone material in mixture, irrespective of 
whether chemical reaction has occurred (—) (equation (3)) 


g mass-transfer conductance for gas phase (Ib. mass/ft.* h.) (equation 
(1)) 


g, constant in Newton’s Second Law of Motion, =4°17 x 10° (Ib. mass 
ft./Ib. force h.?) (equation (16)) 


h_ specific enthalpy of mixture, taking due account of heats of formation 
and reaction (B.t.u./Ib. mass), (equation (2)) 


Hw heat of ablation (B.t.u./Ib. mass), (equation (18)) 
H. heat of ablation without radiation (B.t.u./lb. mass), (equation (24)) 
1,,1, quadrature expressions (—) (equations (47), (48)) 


J mechanical equivalent of heat, Joule’s constant, =778 (ft. lb. force/ 
B.t.u.) (equation (16)) 
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dimensionless parameter measuring relative magnitude of ¢ and of 
the proportionality constant of the phase-changing reaction (Ref. 10) 
(—) (Fig. 3) 


mass transfer rate into gas phase (Ib. mass/ft.* h.), (equation (1)) 


mass transfer rate from solid to liquid phase (lb. mass/ft.’ h.), 
(equation (14)) 


total ablation rate, normally equal to mm’, (lb. mass/ft.* h.), (equation 
(9)) 


mass rate of liquid flow along surface, per unit surface area (Ib. mass 
ft.* h.), (equation (8)) 


molecular weight (Ib. mass/Ib. mole), (equation (28)) 
a number (—) (equation (28)) 


a number indicative of the steepness of the temperature dependence 
of liquid viscosity (—-) (equation (56)) 


pressure (Ib. force/ft.*), (equation (32)) 


heat flow rate per unit area across control surface indicated by suffix 
(B.t.u./ft.* h.), (equation (4)) 


radiative heat flux from condensed phase, supposed to emanate from 
region on side of L surface away from gas (B.t.u./ft.* h.), (equation 
(8)) 


absolute temperature (°F abs.), (equation (17)) 
velocity of liquid interface in x-direction (ft./h.), (equation (32)) 


velocity of gas stream in x-direction, just outside the boundary layer 
(ft./h.), (equation (32)) 


velocity of gas approaching the shock wave (ft./h.), (equation (16)) 
velocity in the y-direction (ft./h.), (equation (38)) 


distance along surface measured from stagnation point (ft.), (equation 
(32)) 


distance normal to interface into the gas phase (ft.), (equation (13)) 
distance measured in y-direction, but measured relative to the solid 
nose-cone material instead of to the interface (ft.), (equation (61)) 

a number (—) (equation (28)) 

a number (—) (equation (27)) 

thermal exchange coefficient (=thermal conductivity divided by 
specific heat at constant pressure), (lb. mass/ft. h.), (equation (38)) 
emissivity of interface (—) (equation (17)) 

dimensionless distance in y-direction (—) (equation (37)) 

time (h.) (equation (13)) 

dynamic viscosity (Ib. mass/ft. h.), (equation (32)) 

density (Ib. mass/ft.*), (equation (13)) 
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€ dimensionless liquid enthalpy (—) (equation (52)) 


Stefan’s radiation constant, =0:1713« 10 [B.t.u./ft.2 h. (°F abs.)*], 
(equation (17)) 


7 shear stress (Ib. force/ft.*), (equation (33)) 


yY transpiration factor (—) (equation (27)) 


Each refers to a mixture state. For explanation, see the text (chiefly 
ST.T. | Section 2.1) and Figs. 1 and 3. 

Ts, 00 


2. The Influence of Thermodynamic Properties on Ablation 
2.1. SOME RELEVANT RESULTS OF MASS-TRANSFER THEORY 
2.1.1. General Remarks 


Figure | illustrates an element of the interface separating a condensed phase 
material from a gaseous phase. A thin control volume, with surfaces § and L, 
encloses the interface element. We shall be concerned with the mass flow rate from 
the condensed phase into the gas, si” (Ib. mass/ft.* h.), and the two heat fluxes g”;, 
and gq”; (B.t.u./ft.? h.); the heat fluxes are unequal, in general, because chemical 
reactions and phase changes occur at the surface. 


Ficure 1. Illustrating notation used in describing heat and mass fluxes at the forward stagna- 
tion point of an axially-symmetrical body when melting ablation may occur. Note that the 
interface between solid and liquid cannot be clearly identified for all materials. 
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Figure 1 also shows, for generality, a second interface separating two regions 
of the condensed phase; two control surfaces, M and N, lie adjacent to it. 
Ordinarily liquid will be present between the surfaces L and M, with solid on the 
other side of N. The control volume of which L and M are opposite faces is 
completed by the surface marked P; material (liquid) is supposed to flow across 
this surface at the rate m”, (lb. mass/ft.? interface area h).* Heat fluxes through 
the control volume are those marked gq”, and g”y, which are conductive in 
character, and g”,,4 which accounts for radiation from the liquid to the gas and its 
surroundings. Conduction through the P surface is supposed negligible but could 
otherwise be added to g” +a... Enthalpy fluxes will be referred to later. 


Also drawn in Fig. | is a control surface O. This is situated within the bulk 
of the solid in a region in which there are no temperature gradients: there is there- 
fore no heat flux through surface O. Conduction through the side-walls of the 
control volume with faces N and O will be neglected. 


It will be supposed that no significant amount of chemical reaction occurs 
within the condensed phases: it follows that the two heat fluxes g”,, and g”,y differ 
only by an amount which corresponds to the latent heat of melting of the material 
crossing the interface between M and N. 


Some materials, for example the glasses, do not exhibit any sudden transition 
from the solid to the gaseous phase. For such materials the surfaces M and N, with 
their corresponding heat fluxes, cannot be identified. 


2.1.2. Gas-Side Conductance and Driving Force 


It is a consequence of mass-transfer theory’: '’:'” that, under certain conditions 
described in detail in Ref. 11, the rate of mass transfer into the gas phase, m’”, is 
equal to the product of a conductance g and a driving force B, namely, 


The conductance is deduced from boundary-layer theory, and the driving force 
from thermodynamic data. Thus the conductance for an axi-symmetrical stagna- 
tion point in steady laminar flow can be read off Fig. 2, which is valid for a Prandtl 
or Schmidt number of 0:7 and for a wide range of temperature ratios across the 
boundary layer’’’. It will be noted that g is a weak function of B. The shape of 
this function must be expected to vary somewhat with the materials in question as 
a consequence of non-uniformities of viscosity, thermal conductivity and diffusion 
coefficient; these influences are, however, of secondary importance to the present 
discussion (see the discussion following equations (27) and (28)). 

In the ordinate of Fig. 2, the symbols have the significance given in the 


notation, with suffix G denoting gases in the main stream just outside the boundary 
layer. 


*In general, liquid will flow out through one part of P and in through another part: this is 
not so, however. in the symmetrical stagnation point flows to be dealt with in the present 
paper 
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nm 


8 
FiGureE 2. Dimensionless mass-transfer conductance D against mass-transfer driving force 
B for the laminar axi-symmetrical stagnation point, with Prandtl and Schmidt numbers equal 
to 0:7, temperature ratios 7,,/T, between 0°25 and 4:0, and property variation laws valid 
for air 


The driving force B in equation (1) is related to the thermodynamic properties 
of the mixtures prevailing in the main gas stream just outside the boundary layer 
(G state), immediately adjacent to the interface (S state), and in the transferred 
substance (T state, defined in Section 2.3.1). Specifically, B is related to the specific 
enthalpy / (B.t.u./lb. mass) and the mass fraction f (lb. mass/Ib. mass defined in 
Section 2.3.1) in these states by the equations 


he —hs 
B h,—hy 
fi—fs 
B= . (3) 
fs—fr 


2.1.3. Definitions 

In equations (2) and (3) the suffixes indicate that the properties must be 
evaluated for mixtures in the states indicated. For the G and S states this presents 
no difficulties. The T state, however, requires special definition. For all the 
problems to be dealt with in the present paper, the following detinition will suffice, 
though it is not the most general : 


fr=l. . (5) 
Equation (4) may be interpreted as meaning that the “enthalpy of the transferred 


substance”, h,, is defined so that, when multiplied by the mass-transfer rate, it is 
equal to the net energy (i.e. enthalpy plus heat) flux across the L or S$ control 
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Ficure 3. Enthalpy-composition (h-f) diagram for a hypothetical nose-cone material mixed 

with air, exhibiting typical air stream (G) and material supply (O) and associated states, 

together with some important heat fluxes. Apart from the ft, (or ¢,) isotherm, all bold lines 
represent boundaries between regions of different phase. 


surfaces. Here suffix 7; means “transferred material in the phase and state of 
chemical aggregation in which it crosses the S control surface.” For generality, 
h, should be replaced by hr,; however, these quantities are equal in all the 
problems considered in the present paper. 


The quantity f stands, in all the examples dealt with in the present paper. for 
the mass of nose-cone material which is required to be mixed with material from 
the gas stream in order to form unit mass of the mixture in question. Thus f=1 
for pure nose-cone material; f=0 for pure gas stream material. Specification of f 
implies nothing about the state of chemical aggregation of the mixture. 


2.1.4. Representation of Thermodynamic Properties of Mixtures on an Enthalpy- 
Composition Diagram 


Another consequence of mass-transfer theory, explained elsewhere’ is that 
the various states already mentioned (G, S, T, T., L. M, N. O, P) can be represented 
as points on a diagram in which the enthalpies, A, of mixtures of the nose-cone 
material with the gas stream material (i.e. air) are plotted against f for various 
temperatures. 


Figure 3 shows such a diagram for a material exhibiting a definite melting 
point, vaporisation or sublimation by surface pyrolysis, and exothermic gas-phase 
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chemical reaction with air. Typical positions of the state-points are shown. The 
following matters are noteworthy. 


(i) State-point G lies on the pure-air line (f=0). 

(ii) State-points L, M. N, O, P, T and T., lie on the pure-nose-cone-material 
line (f= 1). 

(iii) State-point $ lies on a curve, the position of which is dictated by the 
value of a dimensionless quantity K,,,. which measures the speed of the surface- 
pyrolysis reaction relative to the gas-side conductance. 

(iv) Some additional points are marked (7’, G’, G,. G,) for reasons to be 
explained. 

(v) State-point P is placed so that /, represents the mass average of all 
material flowing out through the P control surface. So P lies ordinarily between 
L and M. 

(vi) State-points G, S and T lie on a straight line; this is a consequence of 
the similarity of the differential equations governing the distributions of A and f in 
the boundary layer, provided that the Lewis numbers are unity. 

(vii) The driving force B may be represented by the ratio of the lengths 
GS and ST, i.e. by reason of equations (2) and (3), 


GS 


It will be noted that, from the geometry of similar triangles appearing on Fig. 3, 
we may also write 


(6) 


GS GS GG GG 

Sr TL 

(viii) State-points G, and T, lie on opposite ends of the gas-phase isotherm 

through S. This is because they represent the states of main stream and transferred- 

substance materials respectively, both in the phase of, at the temperature of, and in 
a State of, chemical aggregation appropriate to the S state. 

(ix) Clearly the problem of determination of the mass-transfer driving force 
reduces, when a thermodynamic property chart is available (and when f,,=0 and 
fr=1, as in all the cases here considered) to the determination of the horizontal 
position of point §. Further, since § must lie on a prescribed curve (that for given 
K,,,). knowledge of the surface temperature f. suffices to locate S and so to fix B. 


(7) 


2.1.5. Representation of Heat Fluxes on Enthalpy-Composition Diagrams 


The definition of A;, equation (4), implies that the vertical distance TL, 
measured in B.t.u./Ib. mass, is equal to the quantity q”,/m’”. Since G, lies on the 
intersection of the straight line LS with the line f =0, substitution for TL by way of 
equations (7) and (1) shows that the vertical height G,G is equal to q”,/g. 


The definition of A; also implies that the vertical distance TT, is equal to the 
quantity g”;/rm”, while the vertical distance G.G is equal to the quantity q”s/g. 
The heat flux g”s is, of course, that flowing through the S control surface in Fig. 1. 
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An important quantity is the difference between the heat flux g”, and the 
radiant heat flux q”:.; for this is the heat flux which penetrates beyond the first 
few layers of molecules into the condensed phase. The points G’, and 7” are there- 
fore drawn in Fig. 3, respectively on the f=0 and f=1 lines, such that the length 
G,G’ is equal to and the length 7’L is equal to raa)/m’”. 
Relations (1) and (7) imply that G’, S and 7” lie on a single straight line. 


The points N and M represent respectively the states of solid and of liquid 
mose-cone material at the melting point. The vertical distance between them thus 
measures the latent heat of melting of the material; so NM is equal, when evaluated 
in B.t.u./Ib. mass, to the quantity (g@”’y—q"’y)/m"», where m”, is the mass flux 
through the control surfaces O, M and N in Fig. 1. 


All the foregoing results, together with some consequential ones, are indicated 
in the margins of Fig. 3. 


2.1.6. Some Important Results Derivable from the Steady-Flow Energy Equation 
If the ablation process has reached a steady state. it is possible to derive some 
valuable relationships between the heat and mass fluxes introduced so far. Thus, 
application of the Steady-Flow Energy Equation to the control volume comprising 
surfaces L, O and P in Fig. 1| yields the result 
(hy —ho) + mp (hp—ho) (8) 
If we now define a state Q such that it represents liquid nose-cone material at an 
enthalpy which is the weighted average of those of L and P, we can write 
equation (8) as 


mM o 


where point Q lies between points P and L on Fig. 3. 

We already have two lengths on the A-f diagram involving the quantity 
(g”; —q’raa). Comparison of these with equation (10) enables us to deduce two very 
important results. They are 
GiG’ 

g OO 
m” OQ 
mo TL 


(11) 


and (12) 


Since m’”, is usually the quantity which it is desired to find, while g can be 
calculated from boundary-layer theory, equation (11) provides the solution to the 


problem if the lengths appearing on the right-hand side are known. Here it is 
usually only the establishment of the location of Q which presents any difficulty. 


246 The Aeronautical Quarterly 


/ 


MELTING ABLATION 


2.1.7. Equivalent Results for Unsteady Ablation 


If the surface temperature, mass fluxes and other quantities are varying with 
time, equations (8)-(12) are not valid. Without considering the transient situation 
in general, we shall derive some equations corresponding to those of the last 
section for the case in which material properties at any instant vary only with 
the distance v normal to the surface; the stagnation-point on a nose-cone falls in 
this class. 


For the control volume in Fig. | closed by the surfaces L, N and P, application 
of the First Law of Thermodynamics yields 


L 
GQ’ =m" (hp — hy) + m” (h, — hy) 4 | hp dy 
dé 
N 
where 4 = time (hours) 
p=local density (Ib. mass/ft.*), 
. . . dy signifies integration over the instantaneous thickness of the 
‘y control volume. 


Application of the Law of Conservation of Matter to the same control volume 
yields the relation 
I 


+m | p dy} 0. (14) 
N 


For the control volume comprising the surfaces N and O, together with suitable 
side-walls across which no heat or material flows, application of the First Law of 


Thermodynamics yields 


L 

q’y =m" (hy — hy) + \ lek dy}. (LS) 
oO 


Here it has been supposed that the surface O moves through the solid material so 
that the control volume always encloses the same mass of material. 


PRELIMINARY EXAMINATION OF THE STEADY ABLATION PROBLEM 


We are now in a position to appreciate how ablation rates can be calculated, at 
least in the case of the steady state. We begin by considering some extreme 
situations. 


2.2.1. Material of High Melting Point or Viscosity 


Let us suppose that the whole of the material to the right of surface L in Fig. 1 
is either solid or so viscous that a negligible amount flows outward between 
surfaces L and M. Then #r”,/m”. is zero and m” is equal to #m”,; examination 
of equation (12) shows that points 7’ and O therefore coincide. 
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Now QO, representing the state of the bulk of the nose-cone material, is 
normally given; so is G, representing the state of the air stream just outside the 
boundary layer, together with sufficient aerodynamic information (i.e. “«, Pa, 
du,,/dx) to enable a ¢(B) function to be derived, for example, by way of Fig. 2. 
h., incidentally, is given by 


he =h,,+ Uy? /(28,J), (16) 


where h,, is the enthalpy of the air far upstream of the shock wave adjacent the 
nose-cone (B.t.u./lb. mass), u., is the velocity of the air in this region relative to the 
missile (ft./h.), g, is the constant in Newton’s Second Law of Motion, =4:17 x 
10° (Ib. mass ft./Ib. force h*), J is the mechanical equivalent of heat, Joule’s 
constant, = 778 (ft. lb. force / B.t.u.). 


Since T’ coincides with O, its position is known. If therefore we could find G’. 
S could be found as the intersection of G’T’ with the S-curve for the value of K,,. 
in question; then the problem would be solved. 


How can G’ be found? It lies below G at a distance equivalent to q’raa/ xX. 
Now q” rai is related to absolute surface temperature T;. If we consider the case in 
which the surface has emissivity « and radiates to surroundings at absolute zero, 
the relation is 


where is Stefan’s constant, =0:1713 x 10-* (B.t.u./ft.*h.@F abs.)*). Under other 
circumstances a radiant flux from the direction of the gas may make a contribution 
to q’raa.: then equation (17) must be appropriately modified. 


It is now easy to see how the steady-state conditions can be established by a 
trial-and-error procedure, for example as follows : 


(i) Assume that G’ is identical with G (i.e. neglect radiation). 
Hence locate S§ as the intersection of GO with the S-curve, and 
evaluate B, ie. GS/SO. Hence, from Fig. 2, find a first 
approximation to g. 

(ii) Using the g and 7, values of (i), evaluate g”,,./2@ and so locate 
G’ (second approximation). Then find the corresponding second 
approximation to the position of § (by joining G’ and O). B 
and g. 

(iii) Repeat the process until sufficient convergence has been 
achieved, i.e. until the modifications to the positions of S and G’ 
are negligible from one cycle to the next. 

(iv) Finally, evaluate G’S/SO, i.e. B, and so g, and m’”. 

The foregoing procedure can be carried out extremely quickly. Subsequent 
discussion is therefore simplified by regarding G’ as being given in the data of the 


problem. 
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Summarising, it may be said that the steady ablation rate for a solid or near- 
solid material is very easy to calculate; knowledge of the G’ and O states leads 
immediately to the value of the driving force, whereupon the mass transfer rate is 
obtained through the appropriate standard solution of the boundary-layer equation. 


2.2.2. Material with a Definite Melting Point and a Low Liquid Viscosity 


Another extreme case permitting simple calculation is that in which the 
material of the nose-cone melts at a definite temperature to form a low-viscosity 
liquid. A low viscosity implies that the liquid film must be very thin and so cannot 
have a large temperature difference across it; this means that f, and fy are prac- 
tically identical; so points 1, M and Q coincide on the A-f diagram. 


It is seen that the lowness of the viscosity fixes the location of Q, which is all 


that is needed to permit S and G, to be fixed and 1”,,/ ¢ to be evaluated by way of 
equation (11). Now g can be evaluated through the g¢(B) relation derived from 
boundary-layer theory (Fig. 2); so, since B is now known, 7”,, can be evaluated. 


This is the desired total ablation rate. Clearly it is easy to determine at the same 
time what proportion of this material enters the gas phase, by using equation (12). 


2.2.3. Comparison of the Two Extreme Cases 


Figure 4 illustrates the enthalpy-composition diagram for mixtures of air with 
a nose-cone material exhibiting a definite melting point. The air and nose-cone 
material supply states (G and O states) are specified. A comparison is now made 
of the ablation rates which could be calculated (7) on the assumption that the liquid 
viscosity is very high, (ii) on the assumption that it is very low. 


G’,. G,., S, and L, represent the important state-points corresponding to 
assumption (i); they are determined in the manner already indicated. The ablation 
rate is ¢, G,., G,’/OL,. 


G,’, G,... S, and L, represent the corresponding state-points for the second 
assumption. This time, the ablation rate must be ¢, G,., G,’/OL., in accordance 
with the foregoing discussion. 


Which ablation rate is the greater? Clearly the second, for the following 
reasons : 

(a) OL, is smaller than OL,, as may be seen from Fig. 4. 

(b) G, .G,’ is larger than G, ,G,’, partly because of the lower surface 
temperature which lowers G, and partly because of the higher 
position of G’. The latter effect is caused by the reduced radiant 
heat flux occasioned by the lower surface temperature, and by 
the increased value of g now to be explained. 

(c) g, is larger than g,, because B, is smaller than B,; the latter effect 
may be understood when it is observed that S, lies to the left of 
S,, leading to a smaller value of the driving force for mass 
transfer from the liquid to the gaseous phase. Inspection of 
Fig. 2 completes the argument. 
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Ficure 4. Illustrating the determination of the interface (S, L) states for two extreme cases 


(1) when liquid flow can be neglected (indicated by subscript 1); 
(ij) when the melting point is definite and the liquid viscosity is low. 


It is concluded that a low viscosity is undesirable: it leads to a total ablation 
rate which may be many times as great as that which would correspond to a high 
viscosity; how many times depends, of course, on the properties of the mixtures 
which are displayed by the A-f chart, on the locations of the G and O points, and on 
the magnitude of the radiant effects: there is no general rule. A low melting 
temperature and a low latent heat of fusion are also undesirable. 


2.2.4. The Cases of Moderate Viscosity and of Indefinite Solid-to-Liquid 
Transition 


In practice, the steady-ablation condition of a nose-cone is likely to lie between 
the two extreme conditions just described; the ablation rate will correspondingly 
have an intermediate value. Often, knowledge of the two limiting values will suffice, 
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for example to show that the material in question is unsuitable and so need not be 
further considered. Sometimes, however, more precise knowledge is desired: in 
this case it is necessary to consider the hydrodynamics of the liquid film, as in 
Section 3. 


Some materials which have attracted the attention of designers, for example 
the glasses, do not exhibit a clear-cut change from the solid to the liquid phases; for 
these it is impossible to decide on a melting temperature. The enthalpy-composition 
diagram for one of these materials in mixtures with air is shown in Fig. 5; the 
material is Pyrex glass at one atmosphere total pressure with equilibrium assumed 
to prevail in the § state; the data on which the diagram is based were obtained from 
the paper by Bethe and Adams? and are as follows 


Specific heat of both liquid and solid at constant pressure: 0:25 
B.t.u./lb. mass °F. 

Latent heat of vaporisation: 4450 B.t.u./Ib. mass 

Molecular weight of gas-phase Pyrex divided by molecular weight of 
air: 0-72. 

Equilibrium vapour pressure of Pyrex: exp (14-5 — 83,500/7;) atm. 


For materials of this kind, it is possible to find the lower limit to the ablation 
rate by way of assumption (i) of Section 2.2.1; the upper bound, corresponding to 
assumption (ii), cannot be determined with certainty because of the absence of a 
clear-cut melting point. A hydrodynamic analysis of the liquid layer is therefore a 
necessity for materials of this kind. 


2.2.5. An Example of Pyrex Ablation 


The importance of determining the liquid surface temperature can be seen from 
the following example. 


Suppose that a Pyrex heat-shield at a bulk temperature of 500°F abs. is exposed 
to air at a stagnation enthalpy of 8,000 B.t.u./lb. mass (corresponding to 
ux, = 20,000 ft./sec. approximately). Neglecting radiation from the surface, deter- 
mine limits to the quantity m’’,/¢ ,_,, Where mo is the total rate of disappearance 
of solid material, while ¢, ,, is the value of the conductance corresponding to zero 
mass transfer (from Fig. 2, 2, ,)~0°957 

The upper bound can be found by the foregoing procedure. We find 
T;= 5490 °F abs., B=1-15, 0-596, and so 


Hit” = 0:596 = 0-685. 


This value is lower than that for cooling by means of water (B=7-05), but larger 
than that for the ablation of graphite (B=0-39) in the same circumstances (see 
Ref. 10, where these examples have been considered). 

At a temperature of 5490°F abs. the viscosity of Pyrex is very low; so the 
assumption that the condensed phase does not flow must be invalid. Suppose as an 
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FiGcure 5. h-f chart for Pyrex-air mixture at 1 atm. (Data from Bethe and Adams:?.) 


alternative that the viscosity falls effectively to zero at 4000°F abs., i.e. that the 
Pyrex acts as though this were its melting point. Then measurement on the A-f 
diagram of Fig. 5 reveals 
= 7:92, 
B ~0, 


and so mM" ol 1°92. 


Thus a reduction of the surface temperature to 4000°F abs. brings about an 
11-6-fold increase in the ablation rate. Clearly it is imperative to find out what the 
effective melting temperature really is before deciding whether Pyrex is a suitable 
heat-shield material or not. It is shown in Section 3 how this can be done. 


2.2.7. The Influences of Vaporisation, Surface Pyrolysis and Gas-Phase Chemical 
Reaction 

Before proceeding further, it should be pointed out that we have not had to use 

different procedures according to whether vaporisation is appreciable or negligible; 
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it has sufficed merely to note that the S state-point must lie on a prescribed curve on 
the h-f plane. However, a digression to examine how vaporisation effects have 
made themselves apparent may be desirable. 


Let us return to the example of the last section and consider other effective 
melting temperatures than 4000°F abs. If these are lower than 4000°F abs., the 
point S lies on the line f=0 and drops with temperature by an amount corresponding 
to the specific heat of air; if temperatures above 4000°F are considered, on the 
other hand, the point S moves away from the left-hand border and so rises with 
temperature more rapidly than if it continued along that border. The effect of this 
more rapid rise is that the quantity m”,/2, ., falls more steeply with temperature 
than would otherwise be the case: the cause of the more rapid rise is the finite 
latent heat of vaporisation of the Pyrex (which causes isotherms to bend downwards 
at the S-line), coupled with the increasing mass fraction of Pyrex which is present 
in the S state. These considerations should suffice to allow the favourable effects 
of the endothermic phase change to be perceived. 


Nothing that has been said implies that the S-state gases have to be in 
thermodynamic equilibrium with the surface; they will not be when a pyrolysis 
reaction occurs at the surface; but this merely alters the position of the S-curve on 
the A-f plane. Since a full discussion of this matter, and of the role of the dimension- 
less parameter K,,,, has been given elsewhere’, in the present paper the position 
of the S-curve will be regarded as known. 


Similar remarks can be made about the influence of gas-phase chemical 
reaction between the air and the material emanating from the nose-cone: no explicit 
account has been taken of this phenomenon; yet the procedures described take 
implicit and accurate account of it. A full discussion can be found in the earlier 
literature (e.g. Ref. 13 or 6); here, we merely remark that the existence of gas-phase 
reaction makes itself felt through the curvature which it effects in the gas-phase 
isotherms on the A-f plane: Fig. 3 exhibits this curvature whereas Fig. 5, for 
example, does not. 


2.3. COMPARISON WITH THE WorK OF ADAMS AND OTHERS 


The foregoing sections simply discuss in detail procedures and considerations 
which have already been indicated briefly by the present author”: '”. Since, how- 
ever, a purpose of the present paper is to effect a junction between two streams of 
work on ablation, namely, that just referred to on the one hand and that of 
Sutton”, Bethe and Adams and Lees“ on the other, it is time to make some 
preliminary connections. The survey paper of Adams‘ will be used as the main 
source of information about the second stream of work. It should be mentioned 
that this survey paper uses many results extracted from unpublished reports of the 
Avco Everett Research Laboratory; the present author has not had access to these. 


In the present section of this paper it is appropriate to consider the treatment 
by Adams and others of the role of the thermodynamic properties in ablation, and 
to compare it with that of the present paper. Some remarks will also be made 
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about the solutions to the gas-phase boundary-layer equations which are used in 
the two streams of work. Particular attention will be paid to establishing the 
equivalence of concepts and noting any differences of motivation or principle. 


2.3.1. The Heat of Ablation 


The property of the ablating material which Adams makes the focus of 
attention is its heat of ablation, H.« (B.t.u./!b. mass). This is defined, in the present 
notation, by 

(he —he 

Adams’s notation. 
The numerator of the expression on the right of equation (18) is the heat flux which 
would pass through the S control surface if the conductance were that for zero mass 
transfer, while the denominator is the total rate of removal of material per unit area. 


The reasons for concentrating on H., are two-fold. Firstly, in some simple 
circumstances H.; is an easily ascertained constant which is independent of operat- 
ing conditions; thus, for a system in which the gas stream temperature is not too low 
and in which water is used as a surface coolant (but is prevented from flowing along 
the surface), H.., can be taken as around 1000 B.t.u./Ib. mass, i.e. as approximately 
the latent heat of vaporisation. The second reason is that it is felt that, in any given 
aerodynamic situation, (A, —h..) and ¢, ,, are easily calculated; knowledge of 
then permits quick determination of the important quantity mm”, by way of 
equation (18). 


One might comment that the utility of the Hx, concept is diminished by two 
factors. The first is that, as will be seen, the heat of ablation of most materials is 
strongly dependent on operating conditions, for example on hy, and cannot be 
regarded even approximately as a property of the material alone. The second factor 
is that the numerator of equation (18) (Adams’s q») clearly cannot be calculated 
without knowledge of /.... which implies knowledge of the surface temperature at 
least; so the material properties enter into both sides of equation (18). In these 
circumstances the advantage seems to lie with the more usual mass-transfer practice 
of using the genuine (i.e. operating-condition-independent) thermodynamic proper- 
ties which appear on hA-f diagrams, and then focusing attention on the values of 
mM” o/ 49 Which arise in particular circumstances. 


2.3.2. Graphical Interpretation of the Heat of Ablation 
The variability of H.., with conditions can best be appreciated by giving this 
quantity a graphical interpretation. First, equation (18) is rewritten in the form 


mom 
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On introduction of equation (12), which implies the existence of a Steady state, and 
of equations (1) and (6), equation (19) reduces to 
OG 2p 
Hut (2 
T'L 20) 
There is no length on Fig. 3 which precisely corresponds to the right-hand side 
of equation (20). However, the graphical significance of H.,, can be perceived from 
consideration of the following extreme cases. 


(i) When the surface temperature is low, point § lies well over to the left, T 

and 7” lie very far down on the right-hand border. It follows that T’T, and 7’Q 

are approximately equal. and, since B is small, v, ../g is approximately unity 
Equation (20) then reduces to 

Hin = OE. . 


The physical significance of this equation is that the ablation rate is approximately 
equal to the rate of heat transfer to the surface divided by the enthalpy increase per 
unit mass of melted material. This is exactly true provided that vaporisation is 
negligible (f;=0) and radiation likewise = 0). 


(ii) When the melting temperature or viscosity are high, points T’ and O 
coincide, as has been seen in Section 2.2.1. In this case equation (20) reduces to 


(22) 
g 
Since 2, ,,/g is always equal to or somewhat greater than unity, its exact value 
depending on the value of B and on the appropriate boundary-layer solution 
(Fig. 2), we see that H.,, is always equal to or greater than TT, when no liquid 
flow occurs. 


(iii) If, in the conditions appropriate to (ii), radiation is negligible, O, T and 
T’ coincide, as has been seen. We can therefore write 


(23) 


Comparison of equations (22) and (23) shows that the effect of radiation is to 
increase H., by the amount TT’g, 2. 


2.3.3. “Heat of Ablation Without Radiation,” Hef 


In connection with the foregoing section it is worth noting that Adams 
defines a quantity H.., with the name indicated by the sub-heading; in the present 
notation this quantity reduces to 


He = Hatt - (24) 
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Once again, there is no length appearing naturally on the A-f diagram which 
precisely corresponds to this quantity. Its magnitude and trend can be perceived, 
however, by noting that the quantity inside the bracket is approximately equal to 
T;<T’ when B is small. 


2.3.4. Melting Ablation 


Under the heading “melting ablation” Adams discusses the situation which we 
have treated in (i) of Section 2.3.2; i.e. melting but no vaporisation or radiative heat 
loss. He here refers to the computations of Lees“? who postulated a material with 
a definite melting point and latent heat of fusion. Lees’s expression for the heat 
of ablation is 

H.«=latent heat of fusion 
+ specific heat of solid (t, — fp) 
+06 x specific heat of liquid x (t; — ty). 


This may be expressed in the present notation as 
Hoa = (hy — ho) + 0-6 (h, — hy). ‘ (26) 


Comparing equation (26) with that previously derived for the same process, 
namely equation (21), we see that Lees’s expression can be expressed in geometrical 
language as “Q lies between M and L on Fig. 3”, (this was known already), “at a 
location such that MQ is 0-6 of ML”, (this is the new result). Interpreted physically, 
equation (26) states that the average temperature of the melted material is slightly 
greater than the arithmetic mean of the temperatures of the inner (M) and outer (L) 
boundaries of the liquid film. Of course the temperature of the interface, f,, still 
has to be determined; Lees did this. as we shall do below, by considering the 
hydrodynamics of the liquid layer. 


2.3.5. Subliming Ablation 


By “subliming ablation,” Adams’ means mass transfer into the gas phase 
without appreciable liquid flow; this is the case considered in Section 2.2.1, i.e. that 
for which equation (22) is valid. Adams first draws attention to what he calls “the 
well-known heat-blocking effect”; in the present terms this refers to the non-linearity 
of the relation between mass-transfer rate m” and driving force B (equation (1)), 
which is shown by the variation of the conductance g with B in Fig. 2. Adams 
introduces a transpiration factor ¥, which may be interpreted as the ratio 
2/2, ,,- Citing unpublished reports, he states that laminar boundary-layer 
analyses have yielded a relation which, in the present terms, runs 


g 
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In equation (28), 


Nis a number lying between 0-67 and 0-72, 
@ is a number lying between 0-25 and 0-4, 
M,_ is the molecular weight of the transferred-substance material on 
entry to the gas phase. 

Taking M;r,= 29, which is the value for air, equation (27) yields a curve (or 
rather a band of curves since N is only fixed within limits) which runs close to that 
of Fig. 2. It can thus be regarded as a fairly good approximation to the exact 
solutions on which Fig. 2 is based. The molecular-weight term is presumably an 
empirical one which has been found to bring together solutions for other cases in 
which the Prandtl and Schmidt numbers have been unequal and variable through 
the boundary layer. Equation (27) can therefore be regarded as a useful compact 
formula summarising approximately the exact solutions of the boundary-layer 
equations which have been worked out so far. It must not, however, be accepted 
uncritically in new situations, for example those in which combustion occurs in the 
boundary layer, any more than should Fig. 2 


Adams then derives a formula for the heat of ablation in terms of 2. It is* 


hy —ho+B(he-h, 


Lest 


(29) 
»0 


This may be shown, by a few lines of algebraic manipulation, to be identical with 
equation (22). 

Of course, the use of equation (22) or (29) requires knowledge of the surface 
temperature for the evaluation of A«., and so on. Adams quotes experimental data 
for Teflon ablation which fit these equations within about 20 per cent if the surface 
temperature is assumed to be 800°F. At a value of A, of about 8,000 B.t.u./Ib 
mass the data correspond to a value of the driving force B of approximately 8. This 
value may be compared with the values 0°596, 7-05 and 0:39 for non-flowing Pytex, 
water, and carbon respectively, which were cited? in Section 2.2.5. It is clear that, 
from the point of view of weight of heat-shield, Teflon is a much poorer material 
than carbon, and poorer than Pyrex would be if it could be prevented from flowing 


It is worth noting that the heat of ablation of Teflon varies three-fold over the 
range of conditions covered by the experiments cited by Adams: from about 


This is Adams's equation (8) which reads H,,={C,T,+h,+8 (AA) } 
Presumably the term corresponding to the enthalpy of the base state has been disregarded, 
being small. Adams provides another expression, his equation (13), for the case in which 
vas-phase combustion occurs; these equations are both incorporated in equation (29) without 
the necessity to introduce specific heats and a heat of combustion, both of which are difficult 
to evaluate when manv chemical reactions occur, and when temperature differences are so 
large that dissociation takes place. 

‘It must be remembered that radiation was neglected in calculating the values for Pyrex and 
carbon. In an example worked out by Spalding''®’, the radiation reduces the carbon driving 
torce from 0°39 to 0-174. 
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1300 B.t.u./Ib. mass at low velocities to about 4000 B.t.u./Ib. mass. This fact lends 
support to the comment made earlier concerning the diminished utility of the H.,, 
concept, and suggests that the quantity B is a more satisfactory one for the compari- 
son of materials; for B is directly connected with the weight of heat-shield which is 
needed, as well as being the dimensionless parameter which is characteristic of the 
conditions in the boundary layer, and which is familiar from other mass-transfer 
studies. 


2.3.6. Combined Melting and Sublimation 


Under this heading, Adams“? considers the general case; his expression for the 
heat of ablation may then be written, in our terms, 


hg —ho+(m" {hr —he+ B (he - 


Hur | {(he — hi al 


(30) 
Further algebraic manipulation reveals that this has a significance identical with 
that of equation (20), as is to be expected. 

The evaluation of equation (30) requires knowledge of the ratio of the 
sublimation mass flux #7” to the total mass flux m”,. This is easy to determine 
directly in cases such as that considered by Adams in which the nose-cone material 
consists of a mixture of a volatile material (plastic) with a non-volatile material 
(glass); for then it can be assumed that the ratio m”/m”> corresponds directly to 
the mixture ratio. 

In general, however, evaluation of (30) can only be carried out through explicit 
consideration of the vapour pressure of the subliming material, information which 
is incorporated in the present method into the location of the 5-line on the h-f 
diagram. Adams presents the results of one set of calculations of this kind, namely 
for quartz; experimental data are cited which agree fairly closely with the 
predictions. 


2.3.7. Discussion 

From comparison of the theories of Adams and others and of the present paper, 
the following comments appear justifiable. 

(a) The general formula (20), for the heat of ablation, contains the information 
conveyed by Adams’s various special formulae and by any extensions of these to 
variable specific heats, gas dissociation. etc. 

(b) The use of equation (20) in conjunction with the A-f diagram, with its 
“built-in” S-line, allows the magnitude of the influence of sublimation or vaporisa- 
tion to be perceived more easily than is possible by the use of equation (30); for the 
latter must be supplemented by further equations and numerical computations. 

(c) Neither (20) nor any of the formulae of Adams and Lees are usable for 
determination of H.., unless further assumptions or calculations are made concern- 
ing the surface condition. If liquid flow is absent, however, equation (20) reduces 
to (22), as has been seen, and so does yield H. directly when the A-f diagram is 
available. Adams’s formulae (29), by contrast, cannot be directly evaluated. When 
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liquid flow is present it is necessary to consider the hydrodynamics of the melting 
layer whenever the liquid viscosity cannot be assumed very small. 

(d) In addition to the previously-mentioned matters in which the Adams and 
present theories are in agreement in substance, Adams also uses certain assumptions 
in particular cases for which it is hard to find justification when they are interpreted 
on the A-f diagram. Thus he uses, as a general result, the statement of Lees that 
“the surface temperature for turbulent heating can be roughly 30 per cent greater 


than for laminar heating”. Now the change from laminar to turbulent conditions 


is primarily influential through an increase of the conductance ¢, 59: Yet we can 
see from examination of Fig. 4 that, in vaporisation without liquid flow, for example, 
ts will not change with g, ., at all, if radiation is absent: when radiation is present, 
the effect of ¢, ,, on ts depends greatly on the properties of the material and on the 
location of the G and O state-points. No general statement of this character can 
be made. 

Again, Adams quotes unpublished work of Hidalgo to show that the change 
from laminar to turbulent conditions increases the ratio rm” /m’, by a factor of 
two. Yet clearly such a change is linked to the corresponding change in ¢t;, and the 
nature of the relation depends on the local shape of the S-curve of the A-f diagram. 
among other things. It would be easy to demonstrate that there exist practically 
interesting conditions under which the factor ranges between 0 and 20, say. 

(e) Apart, however, from minor points, such as those just mentioned, 
it can be said that the ablation theory which is advanced by the present author 
does not differ in substance from those of Adams and others, being based on 
the same laws of thermodynamics; it does however differ appreciably in formula- 
tion. The present theory uses graphical presentation of thermodynamic properties 
rather than analytical ones, focuses attention on m’”,/2, ,, tather than the heat of 
ablation, and derives from the practices of chemical engineering rather than being 
invented for a special purpose. 


3. Steady Melting Ablation 
3.1. INTRODUCTORY REMARKS 

We now turn to the problem first considered by Sutton’: the role of 
heat conduction and viscous action in controlling the thickness of the liquid 
layer and so the surface temperature. Sutton wrote down the differential equations 
and solved them numerically for a few cases in which vaporisation was absent. 
Bethe and Adams and Lees”, by making some mathematical approximations, 
obtained more easily applicable formulae, the former for glassy materials, the latter 
for materials with a defined melting point. Other contributors to the theory who 
should be mentioned are Roberts’, Carrier and Lew and Fanucci. The last two 
contributions have not been published, but reference may be found to them in 
Sutton’, 


The following presentation of the theory, although it follows none of the 
authors cited in detail, is guided by their pioneer work, particularly in respect of 
what terms in the equation it is safe to neglect (e.g. liquid-flow inertia terms). 
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Ficure 6. Velocity profile in the liquid film and gas boundary layer. 


3.1.1. Outline of the Theory 

Before plunging into the mathematics, it is as well to attempt a verbal 
description of the processes whose interplay governs the liquid thickness and to 
outline the main features of the theory. 

The thickness of the film of liquid depends hydrodynamically upon the rate of 
melting and upon the surface-shear and pressure-gradient forces which are imposed. 
Thus, for a given melting rate, a large thickness involves a low lateral liquid velocity; 
such a film would be accelerated by the imposed forces. A small thickness for the 
same rate would involve a high liquid velocity and a steep velocity gradient; these 
could not be sustained because of drag at the solid surface. The film thickness must 
therefore have an intermediate value. 

The rate of melting depends thermally on the thickness of the film. For a 
thick film acts as an insulating layer, while a thin one permits heat to be conducted 
across it with only a small temperature difference. 


To the mutual interaction between film thickness and melting rate thus 
revealed, which has many analogues elsewhere in heat-transfer theory (see, for 
examp’*, Nusselt’s™> analysis of the condensation of steam), there is added a 
relation between the heat transfer at the surface and the wall shear. Bethe and 
Adams‘ and Lees‘? employ the Chilton-Colburn’ modification of Reynolds 
analogy here; it will be shown that a further modification is required. This relation 
completes the set of equations needed for a complete solution of the problem; it also 
ensures that the solution is in a form such that the surface temperature depends only 
on the thermodynamic, chemical-kinetic and transport properties of the material and 
gas stream and not on the magnitude of the surface conductance g. 

As a consequence of the fact just mentioned, the outcome of the thermal and 
hydrodynamic analysis of the liquid film is a relation involving the same variables 
as appear in equation (11); simultaneous solution of the two relations, involving 
elimination of the unknown liquid-film state O, leads to the required determination 
of 
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3.2. THE HYDRODYNAMIC EQUATIONS FOR THE LiguID LAYER 


We consider the flow conditions in the liquid film for which the velocity 
distribution is illustrated in Fig. 6. Attention is restricted to the steady state: the 
non-dimensional velocity u/u, is supposed to be a function of the normal distance, 
y, alone. For generality, surface M is shown at the plane at which the liquid 


velocity falls to zero; this plane cannot be identified for the glass-like materials 


volume of radius x formed by the L. M and P contro! surfaces (Fig. 1) leads to 


Mass Conservation. Application of this principle to the cylindrical control 


Momentum Conservation. Application of Newton’s Second Law of Motion to 
a small constant-x element of the control volume in question leads, if the convective 
contributions (!nertia terms) are neglected, to 


d d(u/u;) dP du 
uy Coy Poll (32) 
a\ a\ ax dx 
where the liquid viscosity (Ib. mass/ft. h.), varies with y, 


P is the local jiquid pressure, assumed dependent on x alone 
(Ib. force / fi.*), 
°,=constant in Newton’s Second Law of Motion 


4:17 x 10° (Ib. mass ft./Ib. force h.*) 


On integration twice with respect to vy, equation (32) yields 


u 3 


Vy 


where 7, equals the shear stress at the interface (Ib. force /ft.*) 


Equation (33) may be inserted in (31), leading finally to 


5 ly d 4 
which may be rewritten in dimensionless form as 

\ ug 2¥ ( | Ax) Pr J 

Where A= gor, /UcV (pot: /dx), a dimensionless shear coefficient, (36) 


and »=(y—y,)V (pete a dimensionless normal distance. (37) 
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It should be noticed that the quantity ¥ (poe du;/dx) has already appeared, 
for example in the ordinate of Fig. 2; it has the dimensions of a mass flux. For a 
Prandtl number of 0:7, it is equal to g, ,, /0°957, as has been seen already. 

The double quadrature on the right-hand side of equation (35) cannot be 
evaluated without knowledge of the relations of « and p to », and of the quantity ny. 
To obtain this knowledge, thermal considerations must be introduced. This is done 
in the next section. 


3.3. HEAT CONDUCTION AND CONVECTION IN THE LIQUID LAYER 
3.3.1. The Enthalpy Profile 
Application of the Steady-Flow Energy Equation to an elementary control 
volume within the condensed phase leads, for the flow under consideration, to 
d _dh 
dy\ dy) dy 


where |’ is the local “thermal exchange coefficient” (thermal conductivity 
divided by specific heat at constant pressure) of the condensed 
phase (Ib. mass/ft. h.), and 


v is the flow velocity in the y-direction (ft./h.). 

We here make an approximation, comparable in unimportance to the neglect of 
the inertia terms in equation (32), by replacing pv by m’’,; clearly this substitution 
is correct for regions within the solid phase and in error at the L-surface. With its 
aid we can integrate (38), obtaining, after insertion of appropriate boundary 
conditions, 

,dh 


[ -m”o(h—ho), : (39) 
dy 


which we shall write in dimensionless terms, by way of equation (37), as 


(40) 
(h— ho) dy V dx) 
Equation (40) can be further integrated, yielding 
(p, du idx) d(h hy) (41) 
Jie 
3.3.2. The Convective Flux 
The enthalpy flux in the liquid film is given by 
XM" p 
(Ap — ho) pu(h—ho)dy. (42) 
2 
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Inserting (33), (36) and (37), we obtain 


(' ) (A | (43) 


"M 
This equation enables the enthalpy of the P state to be determined from the 
viscosity and enthalpy distributions. Once h, has been determined in this way. it 
will be possible to determine /,, by using equations (8) and (9), which yield 


m”, 


3.4. RELATION BETWEEN WALL SHEAR AND MASS-TRANSFER CONDUCTANCI 
3.4.1. General Remarks 


The quantity A defined by equation (36) may be determined by examination of 
the solutions of the equation of motion in a laminar boundary layer. Fig. 7 shows 
A obtained in this way as a function of m’”/ /( du,,/dx). The curve is valid for 
the axi-symmetrical stagnation point and is based on a compilation of boundary- 
layer solutions''*?. 


Also plotted in Fig. 7, with the same abscissa, is the quantity D valid for a 
Prandtl number of 0-7; D is defined by 


D=e/ du./dx).. (45) 

So the D-curve of Fig. 7 is a cross-plot of the curve in Fig 
Evidently, since A and D fall off in roughly the same way with increasing mass 
transfer, the ratio A/D should be approximately constant. This is exhibited in 


Fig. 7 by presentation of a plot of A/D. which varies between 1-4 and 2-0 in the 
range of practical interest. 


3.4.2. The “Analogies” 


According to the Reynolds analogy, the ratio A/D should equal unity 
According to the Chilton-Colburn™ modification of this analogy, A/D should equal 
the Prandtl number («/I") to the two-thirds power, i.e. to 0-79 in the present case. 
. either of these values is close to the actual value of A/D, as has been seen. 


The reason for the “failure” of the analogies is that, in the accelerated flow of 
the stagnation point, the equations of momentum and enthalpy transfer are not 
similar in form: the former contains a pressure-gradient term to which there is no 
analogue in the latter. This situation is in strong contrast to the flat-plate flow, for 
example, for which the Chilton-Colburn analogy is well known to be successful. 


Bethe and Adams and Lees’ did not refer to exact solutions of the equations 
of motion for their values of A. Instead they deduced them from values of D, by 
use of the Chilton-Colburn analogy. Their values of A are therefore around one- 
half of what they should be. It should be mentioned, in addition, that the values of 
D are also only approximate in the case of Lees’s work; the D-values of Bethe and 
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FIGURE 7. Dimensionless wall shear (A) and mass-transfer conductance (D) for gas-phase 


laminar boundary layer at axi-symmetrical stagnation point; Prandtl or Schmidt number =0°7 
fluid properties uniform 


Adams, however, agree fairly closely with those used in the present paper. It will 
appear finally that the ablation rate never depends on A in accordance with a power 
in excess of the one-third; the inaccurate value assumed by Bethe and others does 
not therefore have a very serious effect on their ablation-rate prediction. 


w 


COMBINATION OF RESULTS 


5.1. Equations Determining the Interface Temperature 


Before combining the results obtained so far, two assumptions are introduced 
which, although inessential, will simplify the mathematical expressions: p and I’ 
will be assumed uniform within the condensed phase, having the values p, and [’, 
respectively. However, will still be regarded as variable. These assumptions 
allow (35) to be simplified, and (41) to be written 

D h—ho\ 
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» is now eliminated between equations (35), (40) and (46). There results 


D la h h h dh dh 
Ny (4 ))(h—h,) 
The quadrature on the right-hand side of (47) can be evaluated whenever «/«, is 


given as a function of A. 


Equation (47) is of high importance, for it will prove to be the main 
determinant of the interface temperature, which is the quantity needed to fix the 
ablation rate. Equation (47) does not entirely suffice, however, for it contains a 
reference to the average liquid (P) state as well as to state L. These states can, 
however, be related by simplifying equation (43). There results 


(u/ (h—ho) — ho) 


From division of (48) by (47), we obtain 


hp—ho 


= (49 


This equation relates the P and the L states as required, and also permits the 
Q state to be evaluated thereafter 


3.5.2. Evaluation of the Quadratures. Case (a): Definite Melting Point, Uniform 
Viscosity 


We now turn to the evaluation of /,, /., and the ratio /,//,. Considering first 
the case in which « is equal to », for h = hy, and equal to infinity for h << hy, we 
obtain 


2D/A T | 
1,=4A (log. (1/u)}2 (14+ <5 log-(I/éw)} (50) 
) | 
/ 1 — 1 —log, Ey) Hog — loge Em 2! (51) 


h, —ho’ 


In order to establish the orders of magnitude of these terms, it is noted that, 
since the liquid surface (L) and melting (M) temperatures can never be far apart, 
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the dimensionless enthalpy difference €, will usually be close to unity. Then the 
terms containing D can be neglected, signifying that the pressure-gradient term in 
the momentum equation is negligible, and that the velocity profile u(y) is approxi- 
mately linear. The terms in A reduce as follows : 


A ((1 — Em)? (1 


2 | (54) 


The last of these three equations signifies that the mean enthalpy of the fluid 
flowing in the liquid film is nearly equal to A,, the enthalpy of the liquid adjacent 
the interface: h, differs from h, by two-thirds of the enthalpy difference between 
the L and M states. This result is nearly equivalent to that of Lees‘, already cited 
as equation (26). 


The most important points to note about /, and /, for the present are as 
follows : 
(i) /, and /, can be evaluated as functions of (Ay — ho)/(h, — Ao), with 
fixed hy and hy. For h, is the only variable; and specification of 
h, fixes m”/g and m’,/g by way of relations derived in Section 
2.2; and A and D are functions of m”/ 9g, i.e. of B. 
(ii) Both J, and /, increase as £y decreases, approximately at first 
in accordance with equations (53) and (54). 


3.5.3. Evaluation of the Quadratures. Case (b). No Definite Melting Point, 
Temperature-Dependent Viscosity 


Materials like the glasses, which at normal temperatures are super-cooled 
liquids rather than solids, exhibit no definite melting point, but instead have 
viscosities which fall from practically infinite values at moderate temperatures to | 
low values at high temperatures. We shall follow Bethe and Adams in considering 
a class of materials which obey the relation 
hy, — ho 


where n is usually a number considerably greater than unity. 

Of course, no real material follows this law exactly but, since equation (56) 
permits simple evaluation of the quadratures, it will be convenient to force any real 
viscosity law into the framework of (56) by choosing n so that 


(57) 


n- 


(h, ho) 
dh 


the right-hand side of (57) being evaluated for the real viscosity law. 
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After insertion of (56) into the expressions for /, and /,. we obtain 


A _2(D (n+1) 
n 


. 5 
(58) 
A Z2(DIA)T, n+2 
mo/2 pe (n+ (59) 


2 (n+ 3)" 
(n+ 1)? 2(D/A)l, 


o/Z ug + 2)? 


[ 2(D/A)U, n+2 


(60) 


To establish the order of magnitude of these quantities, it is noted that m might 
be approximately equal to 10. Then in most practical cases the second term in the 
square brackets (i.e. that expressing the effect of pressure gradients) will be small; 
since A is around unity, both /, and /, will be around 0:01. Two further points 
should be noted: 7, and /, do not go to zero at a definite value of h,, as was the 
case in the previous example; the ratio /,//, is closely equal to unity, signifying 
that the average liquid-film temperature is not much below that of the interface. 


3.6. PROCEDURE FOR DETERMINING THE STEADY ABLATION RATE 


We are now in a position to lay down a single general procedure for predicting 
the rate of steady ablation of a material which melts, sublimes or pyrolyses, and 
reacts chemically with the surrounding gas. As will be seen, it rests on the mass- 
transfer procedures of Section 2, coupled with the results of Section 3.5. 


Step (i). An enthalpy-composition chart must be constructed for the materials 
and pressure in question. This automatically takes into account sublimation and 
vaporisation, chemical reaction. variable specific heats, dissociation, and so on. If 
surface pyrolysis takes place, the constants of this reaction influence the location 
of the S-curve. 


Step (ii). The G and O states appropriate to the particular problem are located 
on this diagram. 


Step (iti). A series of points along the S-curve are now considered: for each 
such S-point the left-hand side of equation (47) is evaluated through equations (11) 
and (12), together with the recognition that 7”, is equal to m”j—m”. (In evalu- 
ating (11) the relation between the P, Q and L states has to be obtained from (49); 
very often it can safely be assumed that these states are identical, at least in a first 
rough calculation). 


Evaluation of the left-hand side of equation (47) results in a curve, the general 
shape of which is sketched in Fig. 8, with h, as abscissa. At low values of h,, the 
left-hand side is very large, both because of the high value of u, and because of the 
large value 7”,/¢ resulting from the vanishing of the length OQ. Ata value of h, 
marked ft, max in Fig. 8, the left-hand side vanishes; this results from the vanishing 
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LEFT HAND SIDE OF 
EQUATION (47) 


RIGHT HAND SIDE OF 
EQUATION (47 


max 


h 
L 
FicureE 8. Graphical solution of equation (47), leading to determination of A, and thereafter 
nr” 9/8, etc. 


of m”, when O and T”’ coincide; it corresponds to the first of the two extreme cases 
illustrated in Fig. 4 and discussed in Section 2.2.3 

Step (iv). The right-hand side of equation (47) is now also evaluated and 
plotted as illustrated in Fig. 8, for various values of f,. This is a much shallower 
curve, but it rises, if at all, as A, increases. Thus when the material has a definite 
melting point, the right-hand side only has positive values for h, > hy, as is 
disclosed by equation (53). 

Step (v). The intersection of the left-hand side and right-hand side curves is 
noted. This gives the value of /, at which all the equations are satisfied, and so 
determines the sought-after location of S, Q, and so on; the corresponding m”,/g 
value is obtained from (11), leading directly to the value of m’, itself. 

Of course, an actual calculation can be completed satisfactorily without going 
through quite such a long-winded procedure, i.e. without plotting the curves of 
Fig. 8. Essentially the problem is to solve equation (47) numerically, using the 
relations incorporated in the h-f diagram and other data; any convenient numerical 
trial-and-error procedure may be used. In practice the variation of viscosity with 
temperature is so steep that there is little difficulty in establishing the surface 
temperature quite quickly in any given case. 

Comparing the method just described with those of Sutton, Bethe and Adams, 
Lees, and others, it might be noted that the present procedure is more general in two 
respects: firstly, a quite general viscosity law is allowed for by the use of equation 
(47); secondly, no changes are made in the procedure to account explicitly for 
vaporisation or such processes as pyrolysis which have not received attention in the 
references just mentioned. 

Detailed differences between the present procedure and those published earlier 
lie in the more careful determination of the quantity A, which turns out to be almost 
twice that previously assumed, and in the confinement of the numerical work to the 
evaluation of two expressions, namely the left- and right-hand sides of equation (47). 
Some of the differences can be regarded as matters of personal preference, but 
not all. 
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4. The Theory of Unsteady Ablation 


4.1. INTRODUCTION 


Most of Section 2, and all of Section 3, relate to ablation processes which are 
invariant with respect to time. Often, however. the gas conditions and velocity vary 
too rapidly for a steady state to be set up: it is therefore important to know how to 
predict the ablation rate under transient conditions 


It is impractical to do more than broach the subject here: discussion will 
therefore be restricted to the two extreme cases encountered in Section 2.2. namely 


(i) the case in which liquid flow is negligible. and (ii) the case in which the material 
has a definite melting point and a very low liquid viscosity 


If an exact theory were to be presented, even these two cases would present 
complexities owing to the fact that unsteady boundary layers have been little studied 
as yet. However. it appears that most practical problems can be treated with 
sufficient accuracy by a quasi-steady analysis, in which it is assumed that the heat 
and mass fluxes through tine gas and liquid boundary layers depend on the 
instantaneous G,. § an:' states, etc.. and on the instantaneous values of p,.. 1. 
du,,/ dx, ete.. in accordaiice with the same laws as are valid for the steady state 
The treatment to be given here is of this nature: it has the effect of ensuring that 
only within the solid phase is there a transient problem to be solved 


It thus appears that the problem reduces to one of unsteady heat conduction 
within a solid, and that the peculiar conditions of the ablation process are expressed 
solely through the boundary conditions. Since such problems can be solved by 


many well-known exact numerical or approximate analytical techniques, we shall 
not discuss techniques in detail 


4.2, UNSTEADY HEAT CONDUCTION WITHIN THE SOLID 
4.2.1. The Differential Equation 
With reference to a co-ordinate system. z. which is fixed relative to the solid 
material, the enthalpy distribution within the solid obeys a differential equation 
ab ez Ar 
Here z is related to v by dz=dyv—m",/p. (62) 


Where 7”,, now varies with time 4 
4.2.2. The Boundary Conditions 

Far away* from the ablating interface. the solid enthalpy is usually prescribed; 
its value is of course h,. One boundary condition is thus 


—0O, h=No (63) 


*We are here dealing with the problem of a heat-shield which is effectively semi-infinite. The 
problem of the heat-shield of finite thickness can be dealt with by appropriate modification 
to the boundary conditions 
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At the other boundary of the solid, which we denote by N, the boundary 
conditions are 
Z=2 (64) 


~ a Ne 
| (65) 
( 


It should be noted that /A,, gy and zy, are all functions of 6. The quantity z,. in 
particular, obeys the differential equation 


dz, (66) 

dt p 
where +”, is equal to the mass-transfer rate through the N surface as well as that 
through the O surface 


4.2.3. Initial Conditions 


In addition to conditions specified at particular values of z, it is also necessary 
to specify the enthalpy distribution at some initial instant, 6-0. Often h=h», will 
be a sufficiently accurate description. 


4.2.4. Solution of the Problem 


The differential equation (61) can be solved by straightforward numerical 
techniques; for example, the equation can be written in finite-difference form, 
permitting enthalpies at successive instants to be computed in turn. We merely 
have to assure ourselves that sufficient information is available at the boundaries. 


Equations (64), (65) and (66) require the closest scrutiny. Normally none of the 
quantities Ay, g”, and m’”,, will be included in the data of the problem. If, however, 
they can be computed from quantities which are given, the problem is soluble. It 
will be seen shortly that this is the case. 


4.3. CASE (i): NEGLIGIBLE LIQUID FLow 
4.3.1. The Heat Flux q’, 


Since there is no liquid film in the present case, the N surface is close to the L 
surface. In accordance with our convention, the space between them is that from 
which radiation escapes. Therefore 


Now it was shown in Section 2.1 that the heat flux on the right of (67) is related 
to the length G,G’ on the enthalpy-composition diagram. Therefore, in this case, 


Specification of a particular transient problem will usually be given in the form 


of the functions /,;(4) and ¢, ,, (9). From this it follows that knowledge of the 
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FiGuRE 9. Sketch of influence of /i, on q’, deducible from /-f diagram tor fixed materia 
with fixed G state gp, ., and g(B) relation 

thermodynamic properties incorporated in the A-f chart will permit the heat fiux 

q’, to be calculated as a function of A, at any instant #. Fig. 9 sketches such a 

curve, showing that the heat flux falls as the surface temperature (and so Ay) rises; 

this relation is usually a strongly non-linear one 


The data of a transient ablation problem can therefore be quickly translated 
with the aid of an A-f chart, into a family of lines such as those on Fig. 9, with time 
6 as parameter: we cannot yet know which Ay prevails at any particular 4, but we 
can determine the relation between g”, and /, 


4.3.2. The Ablation Rate 


In the present case. 
GS 
m”,=m” — 
ST 


which can be evaluated once G. ¢, |, and Ay are prescribed. Thus, corresponding 
to the g”, (hy. 9) function just referred to, the data of the problem also imply a 
m’., (hy, @) function. Fig. 10 shows a typical curve of #”,, against A, for a 
particular 4; it will be noted that 7”, increases as Ay increases in this case. The 
shape of the curve is of course entirely dictated by the shape of the /-f chart and by 


the relation between g/g, |, and B 


(69) 


4.3.3. Conclusion 


It is now seen that the requirements for the solution of equation (61) have been 
met: at each instant it is possible to evaluate the values of (¢h/¢z), and of dz, /d4 
in terms of the instantaneous value of Ay; the procedures of numerical integration 
can therefore be put into practice without more ado 
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m 
0 
e 
N 
Fictre 10. Typical relation between ablation rate and solid phase enthalpy at the interface at 
ne instant in transient ablation of a material which does not melt. fA, and y [/, ¢,. (du,./dx)] 


are supposed known for this instant 


4.4. CASE (ii): DEFINITE MELTING Low Liquid VISCOSITY 
4.4.1. The Heat Flux q’’, 


In the present case, the knowledge that the liquid film is thin may be applied 
to equation (13); there results 


(Ap — hy) + (hy, — hy) (70) 


+ (hy Ay) =(h h, ye~ (hy 


Since the viscosity is small, states L, P and M coincide, and Ay is prescribed; 
the terms on the right of (71) can therefore be evaluated from the data of the 
problem as a function of 4 alone. Further, a linear relation obtains between 
and 


44.2. The Boundary Conditions of the Conduction Problem 


It thus appears that, in contrast to the situation in Case (7), fy is here 
independent of time and equal to the enthalpy of the solid in equilibrium with the 
liquid, while the heat and mass fluxes are linked together. 


Despite the different nature of the boundary conditions, they suffice to dictate 
an unequivocal numerical solution procedure; the conduction problem becomes one 
with a fixed wall enthalpy but a moving wall, the rate of wall movement being 
connected to the wall enthalpy gradient. Readers familiar with the numerical 
solution of transient heat conduction problems will recognise this one as presenting 
no special difficulty. 
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5. Conclusions 

(a) Steady ablation rates can be calculated straightforwardly usine the 
standard techniques of mass-transfer practice in two extreme cases: (i) when no 
appreciable amount of liquid flow occurs. (ji) when the material has a definite 


melting point and a low viscosity 


(b) When the liquid viscosity has intermediate values, an additional step is 
necessary: this is the determination of the root of the single non-linear algebraic 
relation (47) in which /,, representing the interface condition, is the main variable 


(c) By use of enthalpy-composition diagrams, the problem of «unsteady 
ablation can be reduced to that of solving a single parabolic differential equation 
with a non-linear wall condition; this condition is obtained directly from the data 
of the problem by reference to the /i-f diagram and to the relevant gas-boundary- 
layer solution. 


(d) None of the procedures has to be varied to take account of such processes 
as equilibrium vaporisation, surface pyrolysis, gas-phase chemical reaction, 
dissociation, or variable specific heats; their effects are all taken into account 
implicitly by the enthalpy-composition diagram. 


(e) The author recommends the quantities or as “material 
properties” by which the performance of a heat-shield material can be judged, in 
preference to the “heat of ablation.” These quantities have values of about 017 
for carbon, but more than ten times as much for some of the materials which have 
been discussed in the literature. The values of recommended quantities, like the 


heat of ablation, normally depend on the operating conditions 


(f) The geometrical formulations which have been presented contain, but are 
more general than, the algebraic ones derived earlier by other authors. However, 
more accurate account has been taken in the present paper of the relation between 
shear stress and mass-transfer conductance in the gas boundary layer. Arbitrary 
variations of viscosity and enthalpy with temperature have been considered 


(eg) The concepts and notation of the standard formulation of convective 
mass-transfer theory appear to be sufficiently flexible and apt to accommodate the 
ablation problem without strain. 
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A Study of the Various Boundary Conditions 
for Electrical Analogue Solutions of the 
Extension and Flexure of Flat Plates 


S. C. REDSHAW, D.Sc. and K. R. RUSHTON, Ph.D 


(Department of Civil Engineering, University of Birmingham) 


SUMMARY: The application of electrical Snalogue methods to the analysis of the 

extension and flexure of flat plates is reviewed and the difficulties encountered 

n the satisfaction of the various boundary conditions are discussed. A new 

nethod for treating certain boundary conditions and the operation of the 

electrical analogue is described. New experimental results for two cases which 

present great analytical difficulty, the flexure of a plate with a free edge and a 
plate supported on columns, ure given 


Introduction 


The extension and flexure of flat elastic plates can be investigated by various 
nethods and a method which has recently proved useful is the relaxation solution 
f the partial differential equations which describe the deformation of the plates 


The computational work involved in obtaining a relaxation solution of the 
appropriate equations is long and tedious and the volume of the work becomes 
prohibitive if a very fine relaxation net has to be used. However, pure resistance 
electrical networks, which are analogous to the finite difference form of the partial 
differential equations, have been shown to represent these problems accurately and 


) 


automatically” 


Both the extension and flexure of plates can be represented in two alternative 
forms; one. a fourth-order equation and the other, two second-order partial 
differential equations; the equations used to represent any particular problem are 
dependent on the boundary conditions’. The network used to solve the second- 
der equations is to be described in another paper 


Redshaw’, Liebmann"*, and Boscher and Malavard* have independently 
suggested similar networks for the solution of the fourth-order equations. With 
networks containing about 100 nodes they each investigated several simple 
problems: the boundary conditions of a simply-supported and clamped edge, and 


Received February 1961 


tv 


fugust 196] 


| 

| 


Ss REDSHAW AND R. RUSHTON 


the extensional problem in which the boundary stresses were specified, were success- 
fully represented. The value of the results from these networks is limited by the 
coarse mesh and, accordingly, a network has been constructed, at the University of 
Birmingham, with a mesh interval of approximately one-quarter of its value in the 
original apparatus. However, with such a large network of 1600 nodes the previous 
methods of satisfying the boundary conditions had to be modified to dea! 
successfully with the large number of boundary points 


New experimenta! methods are described in the present paper. Two further 
boundary conditions have been successfully represented. the free edge and that 
appropriate to a plate supported on columns. These last two conditions are ones 
which cannot conveniently be represented by the relaxation technique 


NOTATION 


v.¥ Cartesian co-ordinates 
n.t directions normal and tangential to boundary 
a.b lengths of sides of plate 
hi thickness of plate. mesh interval 
deflection of plate 
E modulus of elasticity 
Poisson’s ratio 
D flexural rigidity of plate 
gq lateral intensity of loading 
Q..Q, shearing forces 
M..M,,M normal moments and twisting moment 
direct stresses 
« stress function 
rop.R.R’ network resistors 
u,v,V,E network voltages 
i network current 
¥ current feeding resistor 
A constant relating physical to electrical quantities 


M.+M 


( or, ) 
Laplacian operator 
V 


biharmonic Operator 


The Aeronautical Quarterly 


| 
} 


ELECTRICAL ANALOGUE 


y 
5 2 
| | 
| 
7 4 8 
| 
FIGURE 1 


2. The Analogy 


Che similarity between the mathematical equations describing the extension 
and the flexure of plates is well known, Richardson’ probably being the first to 
report this. The extension of plates is described by the biharmonic equation 


where ¢ is the Airy stress function, and the flexure of plates by the Lagrange form 
of the biharmonic equation, 


Vw q (2) 
D 
where w is the deflection of the plate, q is the lateral intensity of loading and D the 
flexural rigidity of the plate. Hence, a method of solving equation (2) is directly 
applicable to equation (1). 


The manner in which the electrical network simulates these equations has been 
fully described by Palmer and Redshaw”. The finite difference form of 
equation (2), 

q h* 


4 
20w,-8 w,+2 wat = w, 
n=1 n=5 n=9 D 


‘ (3) 
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E. UPPER NET 
| EX 200 VOLTS 
Pp 
(Moet My) 3 0 VOLTS 
OR 
oy) 
LOWER NET 


REPRESENTS 


V, x 
w OR p 9 


IOOMILLIVOLTS 


FiGuRE 2. Resistance network. 


(where the nodal points are as denoted in Fig. 1), and the equation developed on 
applying Kirchhoff’s Laws to the electrical network defined in Fig. 2 

4 8 12 R’ r R’ 


are identical if the terms in the square brackets of equation (4) are negligible. The 
relations between the physical and electrical quantities are listed in Table I; in each 
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TABLE I 
| Flexure Extension 
,| 
Electrical Quantity | | | 
| elationship Relationship 
Lower net w v=kh v=ke 
{ M=V?2w 
Rh? 
Middle net V _M,+M, V=k, ( v=k 
q RR’h*\ q 
Upper net E=k | 
pp D 1 rp D 


case the relationship between electrical voftage and the required function is 
expressed in terms of a constant k, or k,,. 


Each of the previous authors considered that the term in square brackets 
becomes negligibly small on selecting the values of the resistances such that R and 
R’ are very much greater than r and p. However, with a network containing a large 
number of nodes, even with the maximum practical difference in voltage, this error is 
not negligible, because the neglect of the term in brackets is equivalent to introducing 
an error of about 2 per cent in the applied load. An alternative is to add to the 
upper net voltage another potential equal to the terms in brackets, which in effect 
cancels out these error terms. In practice it is necessary, in the first place, to 
connect the network without a correction, then, from a quick measurement, the 
additional voltage can be calculated and applied. Since this procedure was followed 
in all experiments, the electrical analogue network was made identical with the 
finite difference equation (3). 

3. Boundary Conditions 
3.1. GENERAL 

The equations of the resistance network have been shown to be identical with 
the finite difference equations which govern the elastic extension or bending of a 
plate but, to analyse any specific problem, conditions must be applied to the 
boundaries of the network. These will express either the boundary stresses in terms 
of the Airy stress function or the form of fixity on the edges of the plate in bending. 

The mathematical forms of these restraints are listed in Table I]. Different 
techniques are required to simulate each of the boundary conditions. 

3.2. SIMPLY-SUPPORTED EDGE 

For a simply-supported edge coincident with either the x- or the y-axis the 
boundary conditions are that the normal and tangential moments are both zero, 
with zero deflection along that boundary, 

Le. M,=0, w=0. 
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TABLE II 
Flexural Problem Extensional Problem 
Boundary Functions | Boundary Boundary Functions Boundary 
Defined Condition Defined Condition 
w Simply-supported 
edge 
V2w (Straight boundaries 
only) 
Direct and shearing 
Clamped edge 
stresses defined 
aw (Straight and curved 
(Straight and curved 
souNGaries) > 
on on boundaries) 
rv Free edge 
cn ct 
(Straight boundaries 
w 
(2—») 0 only) 
..3 
On 


Note: n and ¢ express directions normal and tangential to the boundary 


A consequence of zero boundary deflection is that the tangential moment M, is 
also zero 


The achievement of the analogous conditions on the electrical network presents 
no difficulty for they are equivalent to the earthing of the boundaries of both the 
middle and lower nets, 


y=), v=0. 


The validity of this method of representing the simply-supported edge condition 
has been demonstrated by the analysis of several problems concerning rectangular 
plates with simply-supported boundaries under uniform and triangular distributions 
of load. The deflections agreed with theoretical solutions reported by Timoshenko 
and Woinowsky-Krieger™ to within | per cent, and a comparison of moments 
showed no difference greater than 1:5 per cent. 


3.3. THE FUNCTION AND ITS NORMAL GRADIENT SPECIFIED 


The boundary condition fixing both the function and its normal gradient 
frequently arises, for it represents both the clamped edge in the flexural problem and 
also the extensional problem, in which the direct and shearing stresses are defined 
on the boundary. These conditions are expressed in terms of the Airy stress 
function and its normal derivative. 


Unlike the case of the simply-supported edge, there is not one boundary 
condition for the middle net and another condition for the lower net; instead, both 
conditions apply to the lower net and the potentials on the middle net are initially 
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unknown. An iterative method is required whereby the two conditions are satisfied 
on the lower net by applying potentials to both the middle and lower nets. 
Several methods have been suggested to attain these conditions: 
(a) The correct function ~, is applied to the boundary of the lower 
net. The potential of the middle net V, is varied until the slope 
(dv /On), becomes equal to the required value. This method was 
used successfully by Palmer and Redshaw. However, since 
automatic measurement to check whether the conditions are 
satisfied is not possible. this method is only suitable when the 
boundary consists of about a dozen nodes. 

(b) The second method involves setting correctly the gradient at the 
boundary of the lower net by feeding in currents, and then 
adjusting the middle net until the correct potentials on the lower 
net are obtained. This method was tried by Palmer and Redshaw 
but they found considerable instability and cross-coupling, so 
that whenever one adjustment was made all the gradients became 
out of balance. 

(c) Boscher'” suggested that an improved value of the middle 
boundary potential can be calculated from the potentials of the 
lower net. 

Both G. Liebmann‘''’ and the present authors have found 
that, instead of convergence, rapid divergence occurs, even when 
a good approximation to the correct boundary potentials is 
applied initially. 


All the previous methods have suffered from the disadvantage that they require 
the numerical differentiation of the electrical quantities. However, by using current 
feeding resistors the slope is applied as a current and only the value of the function 
at the boundary of the lower net need be measured. 


Both Liebmann and Boscher discussed this method, although they do not appear 
to have fully appreciated its mechanism. 


Figure 3 represents a typical boundary nod From Taylor’s series for the 
nodes 2, 3 and 4, 


ow | (Cow 
But. from Kirchhoff’s Law, 
r 
2 RY y Va), (6) 


where X is the resistance of the current feeding resistor. If the terms of order h 
are neglected the two equations are identical and, in particular, 


ow (7) 
h K, => (u 
OX Jo 
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Considering equation (7), the desired values of 7, and (Ow/¢0x), are known, and 
the value of X specified. Therefore the value of u, can be calculated. This 
potential u, is applied to the current feeding resistor X and, when the correct value 
of v, is obtained by varying the middle net potentials, then both conditions are 
automatically satisfied. 

From the development of this method certain characteristics should be noted. 


(i) The convergence is superior to any of the other methods because, instead 
of one condition being initially enforced, both conditions are made to converge. 

(ii) The sensitivity of this method is critically dependent on the value of the 
current feeding resistor X but, by a suitable choice of resistance, the optimum 
boundary potentials are readily obtained. With a mesh resistance of 100 ohms and 
inter-tier resistors of | megohm the most suitable value for the current feeding 
resistor is 1,000 ohms. 

(iii) Further, this method is equally suitable for curved boundaries. Resistors, 
twice the value of the normal net resistances are required around the boundary and, 
provided that the mesh is sufficiently fine, it is adequate to assume that the boundary 
passes through the nearest nodes. The method of representation is described in 
Appendix I. 

(iv) The fourth desirable characteristic is that it must be possible to check 
rapidly whether the boundary conditions are satisfied. This presents no difficulty, 
since only the boundary potential of the lower net need be examined, and thus a 
method of automatic measurement is feasible. The experimental value of the 
boundary potential on the lower net is compared with the theoretical value set on a 
potential divider, the difference being displayed on a 5 mV. potentiometric pen 
recorder. A uniselector samples each boundary node in turn, and the trace on the 
recorder chart corresponds to the difference between the desired and existing 
boundary potentials. Figure 5, the circuit diagram for a problem discussed in a 
later section, illustrates the method of connecting the recorder to the network. With 
such a rapid means of automatic measurement it is necessary to have a method of 
altering quickly the potential distribution on the boundary of the middle net. It is 
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Ficure 3. Boundary node for simply-supported edge 
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Ficure 4. Simulation of boundary condition of a free edge. 


a satisfactory initial approximation to adjust the middle net potential at every 
eighth node, assuming that the change of potential is linear between these nodes. 
Hence initially, by only altering a few nodes, the potential distribution can be 
changed rapidly; then, when the boundary condition is satisfied to a good 
approximation, the intermediate nodes can be smoothed to fit the general form of 
the curve of the potential distribution along the boundary. Several particular 
problems have been solved satisfactorily. They include the extension of a plate 
subjected to a parabolic stress distribution, the bending of a uniformly loaded plate 
with three edges simply-supported and the fourth built-in, and the bending of a 
uniformly loaded clamped circular plate. A comparison with the theoretical values 
of moment and deflection showed a difference of not more than | per cent. A 
further series of problems entailing this boundary condition is the extension of 
plates containing holes and cracks; these problems have been discussed in another 
paper by Redshaw and Rushton. In the typical problem, discussed in a later 
section, one of the edges was built-in. 


3.4. THe FREE EDGE 


Of the three boundary conditions, that of the free edge is the most difficult to 
represent on the biharmonic network. However, an auxiliary network (Fig. 4), 
has been devised which will simulate the action of a free edge. The theory of the 
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auxiliary network is described in Appendix II. Reference to | ig. 4, will illustrate 
the method of operation. 


(i) Trial voltages are applied to the nodes V., V.. V, . of the 
auxiliary network. 


(ii) The lower net boundary potentials ~. . are set to be 


identical to the network potentials thus 
satisfying the first of the free edge equations 
0x? ey? 
(iii) For the second condition 
+(2-—v) —- . (9) 
CX CXxXCV- 
to be satisfied, the potentials 2,’, . ., Must be equal to 
network potentials 7,,7,.7';.. . . ., Which are one row inside the 


boundary; this check can rapidly be made using the recorder. 
If the potentials are not found to be identical a new set of trial potentials, as (i), 
are applied and steps (ii) and (iii) are repeated. With care, convergence is 
reasonably rapid. 


The simple problem of a uniformly loaded square plate with three edges simply- 
supported and the fourth free was used to demonstrate that this method is feasible; 
agreement hetween experimental and theoretical values being better than 1 per cent. 


If the free edge is adjacent to a simply-supported edge, the auxiliary network 
is a Satisfactory method of simulating a free edge but, with a free edge adjacent to a 
built-in edge, the rapid change in moment at the corner upsets the operation of the 
auxiliary network along the whole of the free edge. For this reason the biharmonic 


network is not suitable for problems concerning cantilever plates. 


4. A Typical Problem 


To demonstrate the mode of operation of the analogue, a problem was con- 
sidered in which each of the three edge conditions occurred. A uniformly loaded 
rectangular plate, with two opposite edges simply-supported, the third built-in and 
the fourth free, was considered. Since the loading was uniform, an axis of symmetry 
occurred midway between the two simply-supported edges. The ratio of the lengths 
of the sides was 2:1, and the shorter side was divided into thirty-two mesh intervals. 
The circuit diagram is to be found in Fig. 5. The experimental procedure was as 
follows: 

(i) An estimate of the moments on the built-in edge was made and 
the corresponding voltage V;, was applied to the boundary of the 
middle net. 

(ii) An estimate was also made of the voltage V; to be applied to 
the auxiliary network. 
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(iii) The voltage on the free boundary of the lower net was set to 
ensure zero normal moment, bank A on the recorder being used 
to speed the operation. 

(iv) The errors for the built-in edge and the free edge were measured 
by banks B and C of the recorder; if all the readings were zero 
the boundary conditions were satisfied. 

(v) Improvements were made to the edge which had the greatest 
error. The middle net potential distribution for this boundary 
was altered until a considerable reduction in errors was achieved. 
Then improvements were made to the other boundary. With 
care this problem can be made to converge in about ten steps. 


The measured potentials which represent the deflection and the sum of the 
moments are recorded in Fig. 6, and the moments calculated by means of the five- 
point formula reported in a previous paper” are plotted in Fig. 7. 


Timoshenko and Woinowsky-Krieger report a theoretical analysis and tabulate 
the moments and deflections at certain points. Good agreement is found between 
the theoretical and experimental results, as shown in Table III. 


TABLE Ill 


Case Theoretical Experimental 


Deflection at the centre of the free edge 0-635 gb4/(Eh3) 0-633 qb4/(Eh3) 
M,, at centre of free edge 0-1172 qb2 0-1196 gb 
Mat centre of built-in edge -0-319 gh2 0-319 gb2 


5. Column Support 


The flat slab supported on columns is of considerable interest in structural 
engineering, for it occurs with both floor and bridge slabs. The respresentation of 
a column by means of the analogue is a simple procedure. 


The effect of a column which is not monolithic with a slab can be considered in 
the following manner. At a certain position in the interior of a loaded plate, an 
upward thrust is exerted until the deflection at that point becomes zero. The effect 
of this upward thrust is identical to that of a column. 


The condition at the point of action of the column can be specified as zero 
deflection with the surface of the plate continuous. An alternative statement of the 
second condition is that the slope of the plate is unaffected by the column. 


In electrical terms, the conditions correspond to zero voltage and no external 
current at the node representing the column. To represent these conditions on the 
electrical analogue a current of opposite sign to the load currents is fed into the 
middle net at the node which represents the column. The current is adjusted until 
the potential on the lower net becomes zero; thus ‘no discontinuity in the deflected 
surface occurs. If there are several columns an iterative technique is necessary. 
The effective shape of the column is a square of side one mesh interval. 
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FIGURE 6 
potentials y and V 


Uniformly loaded rectangular plate having mixed edge conditions 
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Figure 7. Bending moments in a uniformly loaded rectangular plate having mixed edge 
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To verify the accuracy of this method the analogue solution for a uniformly 
loaded square plate, supported on columns and surrounded by identical plates under 
the same loading, was compared with the theoretical solution of Nadai“”. The 
comparison of results is given in Table IV. 


TABLE IV 
Quantity Theoretical ( Nadai) Experimental 
Moment at column 0-245 ga 0-251 ga2 
Deflection at centre of slab 0-0634 gat/( Eh3) 0-0640 ga4/( Eh) 


A further problem which was successfully investigated was that of a uniformly 
loaded rectangular plate simply-supported around the edges and reinforced with 
internal columns. 


6. Conclusions 


The methods described in this paper for treating boundary conditions enable 
many problems concerning the flexure and extension of plates to be solved by 
passive electrical resistance networks. Certain problems, such as that which occurs 
when a free edge is adjacent to a built-in edge, present considerable difficulty. For 
this reason the biharmonic network is not considered to be suitable for the solution 
of problems concerning cantilever plates and an alternative approach is to be 
described in another paper. 


Appendix I 


[TREATMENT OF CURVED BOUNDARIES WHEN THE FUNCTION AND ITS 
NORMAL GRADIENT ARE SPECIFIED 


Provided that the mesh is sufficiently fine, it is satisfactory to assume that the 
boundary passes through the nearest mode. The method of simulation is identical with 
the case when the boundary coincides with either the x- or y- axis, and the method of 
determining the current to be fed into the boundary is illustrated by Fig. 8. 


Appendix II 


THE OPERATION OF THE AUXILIARY NETWORK TO SIMULATE THE FREE-EDGE CONDITIONS 


The two sufficient conditions for the free edge are 


M =0 . (10) 
Cx? oy? 
M..=0. . ‘ (11) 
z ey ry 
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3 4 


4u,+(r/R)V,+2ri,=0 
w, t2w, —4w, —h?V2w, +2h (Ow/ dx), =0 


h(w/)x), =ri, 


(c) Exterior Angle 


BOUNDARY _V, 


20, +20, —4u, +(r/R)V, +4r (i, =0 
+2h /dx+0w/ dy), =0 


5 
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h (Ow / Ox), =2ri,, h(dw/ dy), =2ri, 


FiGure 8. (¢ 


(d) Interior Angle 


2v,+2u +u,-6u,+ 
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2w,+2w +w,—6w,—- 
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h (\w/dx), =2ri,, h(Ow/ dy), 


urved boundaries with the function and its normal gradient specified. 
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First Condition 


Equation (10) can be rewritten as 


dy? 
or, in finite difference form, 
h?M (w,—2w,+w,). . : ; 
In terms of network relationships, this is 
r 
—27 to, ‘ (14) 


If Kirchhoff’s Law is applied to the node 0 of the electrical network (Fig. 4) an 
equation identical with equation (14) is obtained. 


Thus, if the potentials applied to the lower net boundary are such that v7. = 
etc., the first condition is satisfied. 


Second Condition 
The second condition, equation (11), can be written as 
. . (15) 
The two terms Q, and 0M,,.,/cy will be transferred separately into terms of the electrical 
quantities and then equated. 
(a) The Term Q, 
The Lagrange form of the biharmonic equation can be written as 


D 


In finite difference form for the node 0 this becomes 


q oM 
M.+2M.+M.—4M. +h? +2h 0. (17) 
1,+2M,+M,—4M,+h? +2h 


When Kirchhoff’s Law is written for the node 0 an analogous equation is obtained, 


. (18) 
2 4 R 
and, in particular, 
Ox ‘ Rh 


ex 


Hence Q,=D 


It is convenient to convert the voltages V, and V, into terms of the lower net 


voltages; thus 
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(hb) The Term 0M,,,/ cy 


oM 
oy 
Writing this in finite difference form, and using the identity 7 =k,w, we have 
oM,, D 


For the second free-edge condition (equation (15)) to be satisfied, the right-hand 
sides of equations (21) and (23) must be equal. This is achieved only if 7,’= 


Us and SO On. 


REFERENCES 


1. SOUTHWELL, R. V. Relaxation Methods in Theoretical Physics. Clarendon Press, Oxford, 
Vol. 1, 1946; Vol. 2, 1956. 


2. Pa-MerR, P. J. and RepsHaw, S. C. Experiments With an Electrical Analogue for the 


Extension and Flexure of Flat Plates. The Aeronautical Quarterly, Vol. VI, p. 13 
February 1955 


3. RepsHaw, S. C. and Rusutron, K. R. An Electrical Analogue Solution for the Stresses 
Near a Crack or Hole in a Flat Plate. Journal of the Mechanics and Physics of Solids, 
Vol. 8, p. 173, 1960. 


4. Fox, L. and SourHweti, R. V. Relaxation Methods Applied to Engineering Problems 
VIIA, Biharmonic Analysis as Applied to the Flexure and Extension of Flat Elastic 
Plates. Phil. Trans. Roy. Soc. A, Vol. 239, p. 419, 1945. 


5S. Repsuaw, S. C. Electrical Analogues for the Solution of Problems Concerning the 
Extension and Flexure of Flat Elastic Plates. A.R.C. 15,335, May 1952. 


6. LieEBMANN, G. A Resistance-Network Analogue Method for Solving Plane Stress 
Problems. Nature, Vol. 172, p. 78, July 1953. 

Boscuer, J. and MaALavarpb, L. Détermination Analogique des Fonctions Biharmoniques 
Bulletin No. 8, Société Francaise des Mécaniciens, 1953. 

8. RICHARDSON, L. F. The Approximate Arithemetical Solution by Finite Differences, With 
an Application to the Stresses in Masonry Dams. Phil. Trans. Roy. Soc. A. Vol. 210, 
1911. 


9. TIMOSHENKO, S. and Woltnowsky-KrieGer, S. Theory 


of Plates and Shells. Second 
Edition, pp. 117, 133 and 208, McGraw-Hill, 1959. 


10. Boscuer, J. Sur |’ Application de la Méthode des Réseaux Electriques au Calcul de ta 
Deformation des Plaques Elastiques. Comptes Rendus, Académie des Sciences, Paris, 
Vol. 238, p. 1189, 1954 


LiEBMANN, G. The Solution of Plane Stress Problems by an Electrical Analogue Method 
British Journal of Applied Physics, Vol. 6, p. 145, 1955. 


2. Napat. A. Elastische Platten. Springer, Berlin, p. 155, 1925 


292 


The Aeronautical Quarterl» 


| 
| | 

| 


The Static Aeroelastic Deformation of 
Slender Configurations* 


ParT I: Some Elementary Concepts of Static Deformation 


G. J. HANCOCK 


(Department of Aeronautical Engineering, Queen Mary College, 
University of London) 


SuMMARY : By consideration of a series of simple models the static aeroelasticity 

of conventional aircraft is first discussed and then extended to the integrated 

configuration combining fuselage, wing and tail. It is shown that, on the basis 

of linear aerodynamics, a maximum speed exists at which control is possible. 

Overall deformations and control forces become large at speeds approaching 

this critical condition. Part Il, to be published, contains calculations on slender 
plate aircraft including non-linear aerodynamics. 


il. Introduction 


The contemporary language of aeroelasticity was essentially summarised by 
Collar in his classic paper “The Expanding Domain of Aeroelasticity” in 1945. 
This paper unified the whole aeroelastic field by relating aeroelastic effects on 
overall aircraft manoeuvre and stability to the interaction of aerodynamic, structural 
and inertial forces. Static aeroelasticity was distinguished from dynamic aero- 
elasticity by situations where the effects of the inertial forces are negligible. The 
sub-division of these general fields into distinct problems such as divergence, control 
reversal, flutter, effects of aeroelastic deformation on dynamic stability, dynamic 
Tesponse, represents a most significant contribution to the understanding of aero- 
elasticity. The basis of the separation of these various problems is due to the 
different order of frequencies associated with each problem. For example, if the 
frequencies associated with flutter are of a higher order than the frequencies 
associated with overall dynamic response, the coupling between the two problems 
will be weak, so that both problems may be considered separately. Collar suggests 
that static aeroelasticity may be regarded as the class of problem in which its 
associated order of frequency is zero. That these various assumptions are justified 


*A summary of part of this work was read at the meeting of the AGARD Stability 
and Control Panel in Brussels, April 1961. 


Received January 1961. 
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is shown by the success of aeroelastic design calculations over the past twenty-five 
years. But it should be realised that this success has been due mainly to the fact 
that the basic aircraft configuration has remained unaltered during this period, with 
the classic layout of wings, fuselage and tail. However, the forthcoming supersonic 
era may introduce an integrated configuration, combining the wings, fuselage and 
tail into a single structural unit, where the previous ideas based on conventional 
aircraft will not be directly applicable. It is therefore desirable to assess, at this 
Stage, the aeroelasticity of integrated configurations by establishing how far our 
conventional ideas may be extended. It might be expected in such a study that 
there will be a beneficial feed-back in our understanding of the aeroelasticity of 
conventional aircraft. 


Before attempting to discuss any particular aeroelastic effect its seems necessary 
to discuss the term “static aeroelasticity” more fully. Two definitions are given in 
the first paragraph; the first refers to the absence of inertial forces, the second is in 
terms of frequency. These are vague and stem essentially from the simple calcula- 
tions which may be performed for so-called specific static aeroelastic effects. Direct 
extension of these ideas soon leads to difficulty. To clarify these concepts one 
should return to the fundamental problem, namely, the problem of an aircraft in 
flight. Such an aircraft will either be in a steady manoeuvre (i.e. trimmed level 
flight, steady circling flight, steady glide, steady climb) or in a controlled transient 
period from one steady state to another (i.e. during the response to control). In 
addition, during a steady manoeuvre the aircraft must be dynamically stable (static 
stability may be regarded as a particular case of dynamic stability). Thus the total 
consideration of an aircraft’s motion involves two distinct type of calculations; the 
first entails the calculation of the equilibrium of the aircraft in a steady manoeuvre, 
while the second involves the calculation of the dynamic stability of that 
equilibrium. Both of these calculations involve aerodynamic, structural and 
inertial forces. In the steady equilibrium problem the forces (aerodynamic, struc- 
tural and inertial) are independent of time. In the dynamic stability problem the 
forces are functions of time. The problems of transient response are essentially 
identical to the dynamic stability problems; only the boundary conditions are 
modified. Static aeroelasticity is now defined as the domain of aeroelasticity which 
is incorporated into the calculations of the equilibrium of an aircraft in a steady 
manoeuvre. Dynamic aeroelasticity is the domain of aeroelasticity which is 
incorporated into the stability calculations; the further reduction of dynamic aero- 
elasticity into aeroelastic effects on the long period oscillation, short period 
oscillation and flutter will depend on whether or not the orders of frequencies in 
these problems are distinct. The traditional concepts of aeroelasticity all fit in with 
these definitions of static and dynamic aeroelasticity except one, namely, divergence. 
Divergence, as the word literally implies. is a dynamic phenomenon and should 
arise only from dynamic aeroelastic considerations. This is discussed more fully in 
Section 3. 


In this paper discussion is limited to a qualitative assessment of static 
aeroelasticity as defined in the previous paragraph, first reiterating the fundamental 
points for the coventional layout and then investigating the extension of these points 
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to the integrated configuration. Since, from a qualitative point of view, one 
manoeuvre case is as representative as any other, discussion will be restricted to the 
elastic aircraft trimmed for level flight. To illustrate these fundamental points a 
series of simple models is formulated which idealise complementary aspects of 
aircraft elasticity. 


NOTATION 


a,,a_ lift-curve slope of wing and tail 
c,¢Cy chord of wing and tail 
ec distance between aerodynamic centre and flexural axis 
EI beam bending stiffness 
K,.,Ky,,Ky, constants defining the lifts (equations (15) and (16)) 
K constant defining load distribution on lifting beam 
Ly, Ly lift acting on wing and tail 


l, distance between wing aerodynamic centre and aircraft centre of 
gravity 


l, distance between tail aerodynamic centre and aircraft centre of 
gravity 


m_ tail moment coeflicient due to elevator deflection 
M, tail moment about tail aerodynamic centre 
wing semi-span 
V relative free-stream velocity 
VY. maximum trim speed or zero elevator effectiveness speed 
W aircraft all-up weight 
| y(x) deflection function for lifting beam 
z aircraft incidence 
@ wing or tail twist 
¢ downwash at tail 
o weight distribution 
pair density 
n elevator angle 
A= —[KpV*/El]'* 


mo torsional stiffness of wing or tail 


2. Model 1 


An attempt is made in this problem to discuss the fundamental ideas of wing 
deformation on a conventional type of aircraft, consisting of wings, fuselage and tail 
unit. It is assumed that the wings, fuselage and tail are all rigid, but that the 
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mounting of the main wings to the fuselage is elastic, allowing the wings to have the 
single degree of freedom of twist only (no bending flexibility). The torsional stiff- 
ness of the wing (mean chord c, span 2s) relative to the fuselage is denoted by 
me (Ib. ft./rad). The angular deformation of the wing relative to the fuselage is 
denoted by @. 


It is assumed that the wing takes up the attitude 4=0 at zero forward speed. 
This configuration is shown in Fig. 1. 


TORSIONAL STIFFNESS ft/rad) 


FUSELAGE 
TAIL 
WING 
Ficure 1. Flexible wing model 
In steady level flight, to a first approximation, 
wing lift-= weight = constant (1) 


This equation implies that the lift on the tail may be considered small compared 
with the lift on the wing in this case. Also the rate of loss of fuel is neglected. The 
aircraft is assumed to be trimmed by the elevator setting of the tail unit at all 
forward speeds in level flight. Thus, from equation (1), 


where a, is the lift-curve slope of the wing. It is implicitly assumed in equation (2) 
that the angle of zero lift of the wing is zero. The structural equilibrium of the 
wing is given by the condition that 


(Lift) (ec) = mA, . (3) 


where ec is the distance of the axis of twist aft of the effective aerodynamic centre, 
assuming that ec > 0. Therefore 


Substituting equation (4) into (2) gives 


W 


kpV*2sc a,me (me — 4pV72s ec*a,) (5) 
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Equation (4) states that the wing angular deformation is constant, independent of 
forward speed. This follows physically from the initial assumption that the lift 
remains constant, and the lift uniquely determines the angular deformation. 
Equation (5) gives the value of 2, the fuselage incidence for level flight at a given 
forward speed V, which is controlled by the elevator at the tail; then z=0 when 


Me 
(6) 
4p2s ec’a, 

However, this speed has no special significance from the deformation point of view 

(4 is constant for all speeds), although it is recognised as the classical wing divergence 
speed. Presumably, for the problem considered so far, increasing the speed above 
that of equation (6) merely results in z becoming negative. Classical divergence 

does not come out of this deformation calculation; it is a separate problem. 


It can be shown that if the present model, with the flexibility concentrated at 
the wing root, is replaced by a more representative model with a flexible wing (i.e. 
flexibility distributed along the span), rigidly fixed at its root chord to the fuselage, 
then the deformation function 6 (y) (where y is the spanwise co-ordinate) does vary 
with forward speed. When the speed is reached for z2=0, the deformation function 
is finite; this is usually quoted as the “divergence speed.” However, further increases 
in speed in this case predict that 2—»-©O at a finite velocity, resulting in 
6 (vy) —> + 0; this differs from the limiting tendencies in equation (5). 


3. On Wing Divergence 


The terms of reference laid down in the Introduction restrict the discussions to 
static aeroelastic problems, but a digression might be permissible to mention a few 
points on wing divergence, which is classified in the literature as a static aeroelastic 
problem, although it is not admissible in this classification under the present 
definitions. This will be illustrated in the following paragraphs. 

The usual argument is that wing divergence is a critical speed above which the 
wing is unstable. By reference to Fig. 1 we consider the incremental moments which 
act on the wing due to a small displacement 64. The incremental aerodynamic 
moment (4pV72s ec*a,56) tends to increase 64, while the structural moment (m6) 
tends to reduce 64. Thus, if 


the initial disturbance 44 will be reduced. But, if 


4pV*2s ec’a, > me, (8) 


any initial disturbance 66 will be amplified. The wing divergence speed is then 
defined in this case as 
Me 


V? (9) 


bp2s ec’a, 
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Equation (9) is identical to equation (6), that is, the speed for 2=0. Now equation 
(9) gives the divergence speed of the wing if the aircraft is restrained not to move 
away from its trimmed condition in level flight during the application of a general 
disturbance 66. This is, of course, an artificial constraint which cannot be imposed 
in practice, where such a general disturbance 66 on the wings will induce a disturb- 
ance of the whole aircraft. Thus equation (9) only has a meaning if the displaced 
wing motion can be separated from the overall aircraft stability motions. But wing 
divergence will necessarily induce overall aircraft divergence, so that coupling 
between wing response and aircraft response must be strong. For a complete 
analysis one should therefore refer to the full longitudinal equations including wing 
flexibility. However, to illustrate qualitatively the foregoing remarks, we shall 
reduce the order of restraint on the aircraft from that implied in equation (9) by 
considering the static stability of the simple model shown in Fig. 1. Suppose that 
this model is pivoted at the centre of gravity, as shown in Fig. 2, and the aircraft is 
displaced through an angle 6z and the wing is displaced through an angle 46 from 


FiGureE 2. Divergence model. 


its equilibrium position. Divergence is defined as the speed at which the whole 
system is in a state of neutral equilibrium, with both 6¢ and 62 indeterminate. The 
neutral equilibrium of the wing is given by 


The neutral equilibrium of the aircraft is given by 
*S wa, (66 + 62) [CL 02/02) 62 (11) 


where ¢ is the downwash angle at the tailplane due to the main wing. When 62 is 
zero, equation (10) gives the classical wing divergence speed. If 66 is zero in 
equation (11), we have the condition that the aircraft with wings rigidly mounted to 
the fuselage has zero static margin. For the present purposes it is assumed that the 
rigid system is statically stable, so that 
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Elimination of 64 and 62 from equations (10) and (11) gives the divergence speed of 
the whole system as 


be (1 S,a,l,e¢ -1 
boSwec a, 1+ : (13) 
(static margin), 
(l S 
c C2 


With typical values of /,/c=0-20 and a 10 per cent static margin, the divergence 
speed is down by a factor of about two compared with the classical wing divergence 
speed given by equation (9). 


Wing divergence could only have some independent meaning if the static 
margin of the rigid aircraft were infinite, which is outside the realm of practical 
configurations. 


4. Model 2 


In Model | the deformation of a flexibly mounted wing on an otherwise rigid 
air frame was considered. Attention is now turned to the “flexibility” of the 
fuselage, to study its deformation characteristics. Again, to illustrate the main 
principles, an idealised version of the problem is taken. The fuselage is represented 


| 
LB/F T 


+ 


Ficure 3. Flexible fuselage model 


by an elastic beam of uniform cross section (i.e. uniform weight distribution 
o Ib./ft.). The lifts on the rigid wings and tail (L. and Ly respectively) act at the 
ends of the fuselage. The part of the fuselage ahead of the wings and aft of the 
tail is neglected. The wing section is assumed to be symmetrical (i.e. angle of zero 
lift is zero) with its centre line parallel to the local fuselage centre line (i.e. the 
incidence of the wing is the same as the incidence of the front of the fuselage). This 
is shown in Fig. 3. 


The overall equilibrium of the “aircraft” is given by 
Ly=L,=el/2. (14) 
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We are effectively neglecting, at this stage, the weight of the wings. This is a 
reasonable assumption in this case since interest is confined to the fuselage 
flexibility. The tailplane is assumed to have the incidence of the fuselage at the 
tail, so the elevator must be arranged to give a lift on the tail equal to the lift on 
the wing. 


If the Cartesian axes (x, y) shown in Fig. 3 are used (i.e. x is the stream 
direction, y is normal to the stream), then we may write 


Ly = —KypV? (dy/dx),-, ERS) 

and L, [ Ky (dy dx), Ky on). ° (16) 

where K,,, Ky,, Ky, ate constants depending on the geometries of the wing and tail, 
and » is the elevator angle. The negative sign is included to counteract the negative 


sign of dy/dx. Therefore, if the fuselage is rigid, 


—dy/dx=a=constant 


ol 
where a= (17) 
ol 1 
and n= Kz,2) (1B) 
For a flexible fuselage with uniform bending stiffness E/, 
Ely”” = 
with boundary conditions 
Neng 
| 
The solution of equation (19), satisfying equation (20), is 
cis #71 
which, on satisfying equations (15) and (14), gives 
= = = 22) 
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The arbitrary constant a is not significant; it implies that the system may be dis- 
placed normal to the stream without changing the essential features of the problem. 
The constant b (=el/(2K,.pV*) does not contribute to the deformation of the 
fuselage, since it represents a pure rotation of the whole system. Thus, from 
equation (22), the deformation of the fuselage is constant and independent of 
forward speed. This is a similar result to that obtained in Model 1. Again, this is 
not unexpected since the loads on the wing and tail are maintained at a constant 
value for all speeds. It follows that the incidence at the front of the flexible fuselage 
(i.e. x=0) is the same as the incidence of the rigid fuselage. The deformation is 
sketched in Fig. 4. 


| FLEXIBLE 
= 
RIGID SS DEFORMATION 


INDEPENDENT 


FiGurE 4. Deformation of flexible fuselage 


We can now return to a point mentioned earlier in this section. If constant 
loads are added to the system at the nose (x =0) and tail (x=/) to simulate the wing 
and tail weights, the main conclusions are unaffected. 


5. Model 3 


To complete the analysis of the static aeroelastic properties of the trimmed 
conventional configuration we must discuss the deformation characteristics of the 
tailplane which is flexible relative to the rest of the airframe. So, for this model, we 
take a rigid wing, fuselage and tail with a flexible mounting of torsional stiffness 
me lb. ft./rad. between tail and fuselage, as shown in Fig. 5. 


The attitude of the fuselage and wing is denoted by z. The angular deforma- 
tion of the tail relative to the fuselage is denoted by 4, where 4=0 corresponds to 
the zero tail force. Now for overall trim, the lift on the tail is maintained at a 
constant value, Ly, say. The moment acting on the tail is usually neglected in the 
overall stability and control calculations, although it is fundamental for the tail- 
plane equilibrium. Now the lift on the tail may be written 


Lr= [a, a.y]=constant. . ‘ (23) 
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FiGure 5. Flexible tail model ‘ 


The moment on the tail about the tail aerodynamic centre M, (positive nose-down) 
is given by 


The structural equilibrium of the tail is given by 

where ec; is the distance of the axis of twist aft of the aerodynamic centre. From 
the overall trimmed condition it must be remembered that 4pV?(2—<) is constant 


Therefore, from equations (23), (24) and (25), the elevator angle to trim is given by 


{me [Lr (2 — — 
{me — 4pV*Srcx (a,m/ a.)} 


(26) 


Therefore » tends to infinity as V tends to V., where 


Me 
(a,m/ a. ) 


This speed V, is therefore the maximum speed theoretically possible for trimmed 
level flight. The elevator effectiveness for trim is defined as the elevator required to 
trim with a rigid tail (/.e. equation (26) with mz infinite) divided by the elevator 
required to trim with a flexible tail. Thus 


Me — / ae) 
4pV*Syecya 
Ly —4pV?S;a, (2-2) 


elevator effectiveness 


me — 


Thus V. may also be regarded as the speed where there is no control effectiveness 
This is sometimes referred to as the elevator reversal speed, but the definition is not 
really admissible, since the aircraft can never attain the speed at which it cannot be 
trimmed and so the cencept of reversal of control, which implies flight through the 
critical region, is not relevant. In practice, however, the maximum speed in level 
flight will depend on the maximum elevator setting as determined by the stalling 
characteristics of the tailplane. 
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The rate at which the effectiveness decreases to zero, according to equation (28), 
depends on the parameter A, where 


ea, L, | 


A 
m 


(29) 


which is independent of the forward speed V. When A=], there is no loss of 
effectiveness for V << V., although at V=V, any elevator angle will trim the 
aircraft. 


6. Model 4 


Model | is concerned with a flexible wing, Model 2 with a flexible fuselage, and 
Model 3 with a flexible tail. An idealised integrated system is now considered in 
which the wing, and fuselage and tail are combined together, as they would be for a 
slender configuration with the aerodynamic loading distributed over the whole of 
the length. 


Ficure 6. Lifting flexible fuselage 


The model is taken to be a simple beam of uniform cross section (uniform 
bending stiffness distribution E/, weight distribution « Ib./ft., total length /) under a 
simplified aerodynamic load distribution which is assumed to be proportional to the 
local incidence. The load distribution is assumed to be KpV~2(x) and K is taken 
to be constant. To simulate the force on the elevator control, a force P is applied 
at the trailing edge of the system, ensuring that the overall equilibrium conditions 
(lift equals weight, and total moment zero) are satisfied. The system is shown in 
Fig. 6. The overall equilibrium conditions are 


dx+P=cl (30) 
and (Kev o)(I-x) dx=0. 
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The deformation function y (x) satisfies the differential equation 


Ely”” = ; ‘ 
The boundary conditions are 
=() 
|. leur =0 


It can be shown that any solution of equation (32) which satisfies boundary 
conditions (33) automatically satisfies the overall equilibrium equations (30) and 
(31). This is physically obvious since the static deformation problem can only be 
solved if there is no dynamic problem due to the out of balance of any resultant 
forces acting on the system. The solution of equation (32) is 


KpV 


y= , + 2 cos 


~ 


where A (KpV*/ED)'* and A, B,C are arbitrary constants. The three boundary 
conditions 


become 
B /3 

0 
B 

+*=C =0 


(35) 
The fourth boundary condition 
y”’=—-P/El 
then gives the force P in terms of A, Band C. Thus the general solution is 
A=B=C=0 
together with 
P=0. ‘ (36) 
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However, according to equation (35) a critical speed exists, when the coefficients 
are not uniquely defined. This critical speed is determined from the determinantal 
equation 


1, 


1, 


3 3 3 

— cos Al + /3 sin Al), —4e-™ | sin $y) /3cos 

(37) 


In general a finite solution of this determinantal equation exists. Thus, at a critical 
speed, equation (34) (for y’) can be reduced to an equation with one arbitrary 
constant, A, say. Then the trailing edge force P can be determined, also in terms 
of A, from the condition 


P= —-Ely””|, 


Thus the general condition for equilibrium is given by 


with P=0, up to a critical speed V., given by equation (37). At V, an alternative 
mode of equilibrium is possible, the amplitude of the mode depending on the 
magnitude of the elevator control force P. 


Now in this particular problem there is no deformation of the system for speeds 
below V. (since y»’=constant). This is due physically to the fact that the 
distribution of aerodynamic loading exactly balances the weight distribution. 


7. Model 5 


Although Model 4 has shown that a critical speed exists at which there are 
infinite positions of equilibrium. the model is restricted because of the unique 
correspondence between the weight and aerodynamic load distributions. For 


ver \ 
— | 
Wo 


FiGure 7. Lifting flexible fuselage. 
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further appreciation of this effect, and at the same time to represent practical 
conditions a little more closely, Model 4 is now reconsidered with a more complex 
weight distribution. This will ensure that the weight distribution and load 
distribution will not exactly balance. 
Model 4 is now reconsidered with a weight distribution given by « {2 —(x/))} 
(instead of the uniform « Ib./ft.). This particular weight distribution satisfies the 
condition of overall static stability if the model is rigid. Again a trailing edge force 
P is included to trim the model for level flight. This is shown in Fig. 7. 


The overall equilibrium is given by 


1 
dy x 
| dx+P= | (2 (38) 
0 
and | | Kev (2 xydx=0. . . (39) 
The differential equation for the equilibrium deformation y (x) of the model shown 
in Fig. 7 is 
Ely’” = —KpV?y' -o (2- 
together with the boundary conditions 
and len = — P/E. (42) 
The solution of equation (40) is 
y=- 1B cos 5 Ax+C sin (43) 


where A= —(KpV?/EI)'* and A, B,C are arbitrary constants. The three boundary 
conditions (41) become 


B /3 
0 
B V3 o 
A- ~ 
2 2 KlpV?a 
- *(cos Al+ 4 3 al) Se 2(sin¥ - /3 cos = 4 


~ KlpV7A 
. (44) 
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The solution of equations (44) gives a set of unique values of A, B and C for any 
given speed V. Then the control force P is calculated from equation (42). 


Now at a critical speed V., given by the determinantal equation (37), the 
coefficients A, B and C, together with the control force P, become infinite. Thus the 
aeroelastic deformation becomes infinite as V—>V., together with an infinite 
elevator force. The large aerodynamic forces (which necessarily give rise to the 
large deformations) cancel the large control force P so that the overall equilibrium 
(i.e. lift equals weight, total moment zero) is still preserved. This loss of control 
effectiveness in level flight, due to the overall elasticity of the structure, is equivalent 
to the loss of effectiveness of the tail of the conventional aircraft described in 
Model 3 (Section 5). 


It is worth noting that the critical speed V. is the same for Models 4 and 5. 
This means that the limiting velocity V. does not depend on the weight distribution, 
but only on the aerodynamic and stiffness loadings. The control effectiveness is, 
of course, dependent on the weight distribution, in addition to the aerodynamic and 
stiffness distributions, and can be defined, relative to the rigid aircraft, in the same 
manner as for the conventional aircraft. 


8. Discussion 


By a consideration of a series of simple and idealised models an attempt has 
been made to spotlight some of the fundamental points of static aeroelasticity of 
both conventional and integrated configurations. Models 1, 2 and 3 together 
illustrate some of the main points concerning the conventional layout in level flight. 
The separate effects of wing (Model 1), fuselage (Model 2) and tail deformations 
(Model 3) may be superimposed for the overall picture. In trimmed level flight 
wing deformation in general varies with forward speed (it does not become infinite 
at a classical divergence speed), fuselage deformation is essentially independent of 
speed, while there is a loss of control effectiveness to trim with increasing speed, 
due entirely to tailplane flexibility. A critical speed is reached at which the aircraft 
cannot be trimmed, which by definition corresponds to the speed of zero elevator 
effectiveness for trim. 


Models 4 and 5 show that an integrated configuration of fuselage, lifting and 
stabilising surfaces suffers from loss of control effectiveness, again with total loss 
of effectiveness at a theoretical maximum trim speed, but this time involving very 
large deformations of the whole airframe. Control effectiveness to trim is related to 
the out of balance between the weight distribution and the aerodynamic loading. 
The speed of zero effectiveness is independent of the weight distribution and involves 
only the aerodynamic and structural stiffness forces. It is possibly desirable to 
calculate the maximum trim speed without reference to the whole panorama of 
aeroelastic deformation against speed. This may be done by reference to Models 4 
and 5, where the maximum trim speed may be defined as the speed at which an 
arbitrary small controlled deformation may be added to the equilibrium state. 
Controlled in this sense implies that this incremental deformation is balanced by 
a small trailing-edge control force so that the total overall equilibrium conditions 
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(lift equals weight, total moment zero) are still maintained. The differential equation 
for a small deformation 6y (x) for the lifting integrated system, as in Model 5, 
is therefore 


with the boundary conditions 
Oy” | = Sy” | . (46) 


where oP is the additional trailing edge control force necessary to maintain overall 
equilibrium. A solution for 6y from equation (45), combined with the three 
boundary conditions (46), will give the critical speed in the usual eigen 
value manner. 


This note has been restricted to the steady manoeuvring case of trimmed level 
flight. The conclusions, however, apply to all cases of steady manoeuvring flight 
(steady circling flight, steady glide, steady climb) when linear theories apply. The 
only difference between a general steady manoeuvre and the trimmed case con- 
sidered here is in the absolute values of the forces normal to the aircraft direction, 
but the distribution of loading will be the same at any given airspeed. Thus the 
speed of total loss of control effectiveness is independent of the manoeuvre. Control 
effectiveness also is a function of speed only, not the manoeuvre (i.e. the normal 
acceleration). 


In this report linear aerodynamics have been assumed. This assumption is 
valid below the stalling condition for a conventional aircraft, which implies that 
the true maximum trimmed speed for this type of aircraft is limited well below the 
theoretical maximum by the stalling conditions over the tailplane-elevator system. 
However, the aerodynamics of the slender integrated configurations are dominated 
by non-linear aerodynamic forces, due to the appearance of leading edge vortices. 
The inclusion of these effects in more detailed calculations of static deformations of 
slender configurations appears in Part II of this paper. 
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